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Abstract 


Numerical simulation of two classes of unsteady flows are obtained via the Navier- 
Stokes equations: a blast- wave/ target interaction problem class and a transonic cavity 
flow problem class. The method developed for the viscous blast- wave/target interac- 
tion problem assumes a laminar, perfect gas implemented in a structured finite- volume 
framework. The approximately factored implicit scheme uses Newton subiterations to 
obtain the spatially and temporally second-order accurate time history of the interac- 
tion of blast-waves with stationary targets. The inviscid flux is evaluated using either 
of two upwind techniques, while the full viscous terms are computed by central differ- 
encing. Comparisons of unsteady numerical, analytical, and experimental results are 
made in two- and three-dimensions for Couette flows, a starting shock-tunnel, and 
a shock-tube blockage study. The results show accurate wave speed resolution and 
nonoscillatory discontinuity capturing of the predominantly inviscid flows. Viscous 
effects were increasingly significant at large post-interaction times. 

While the blast-wave/target interaction problem benefits from high-resolution 
methods applied to the Euler terms, the transonic cavity flow problem requires the 
use of an efficient scheme implemented in a geometrically flexible overset mesh envi- 
ronment. Hence, the Reynolds averaged Navier-Stokes equations implemented in a 
diagonal form are applied to the cavity flow class of problems. Comparisons between 
numerical and experimental results are made in two-dimensions for free shear layers 
and both rectangular and quieted cavities, and in three-dimensions for Stratospheric 
Observatory For Infrared Astronomy (SOFIA) geometries. The acoustic behavior of 
the rectangular and three-dimensional cavity flows compare well with experiment in 
terms of frequency, magnitude, and quieting trends. However, there is a more rapid 
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decrease in computed acoustic energy with frequency than observed experimentally 
owing to numerical dissipation. In addition, optical phase distortion due to the time- 
varying density field is modelled using geometrical constructs. The computed optical 
distortion trends compare with the experimentally inferred result, but underpredicts 
the fluctuating phase difference magnitude. 
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Chapter 1 
Introduction 


Performance envelopes of vehicles and structures which interact with the atmosphere 
are often limited by unsteady aerodynamic effects. Specific examples which are ad- 
dressed here range from the blast- wave/ target interaction problem, where peak over- 
pressures are many times quiescent conditions, to the seeing problem, where density 
fluctuations contribute to image degradation. Only recently have advances in com- 
puting power and numerical algorithms provided the potential, complementary to 
experimental studies, for the timely design of effective configurations. However, the 
use of unsteady computations in the design phase is presently at an immature stage 
of development. The objective of this effort is to demonstrate computational tech- 
nologies as applied to current topics of interest in the unsteady, compressible perfect 
gas regime. It is hoped that, through comparison to accepted experimental data, 
computational methods can make significant contributions to the design of systems 
which interact with unsteady flows. 

In these studies, two solution methods to the Navier-Stokes equations are pre- 
sented and applied to several test cases. The cases pertain to the blast- wave/ target 
interaction or the transonic cavity flow problem classes. The method used to model 
the strongly unsteady blast-wave flows concentrates on resolution of the complex 
physics through the use of characteristic-based schemes. In contrast, for the tran- 
sonic cavity flow problem, the combination of complex geometries and large problem 
size requires the use of an efficient integration scheme. Comparison of the numerical 
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and experimental results will provide a reference base of unsteady numerical results 
for use in configuration design. 

1.1.1 Blast- Wave Problem 

The study of the effects of blast-wave impingement upon vehicles and structures is 
of practical consideration in the determination of their survivability. The experimen- 
tal study of the blast- wave/target interaction problem requires the use of expensive 
above ground tests or facilities such as the U.S. Army Large Blast /Thermal Simulator 
(LB/TS) facility depicted in Fig. 1. Moreover, experiments can suffer from limited 



Figure 1: Proposed Large Blast/Thermal Simulator facility [1] 


phase durations, and deduction of the physics of the flowfield is difficult because 
of practicalities in data acquisition methods. Visualization of the propagating wave- 
fronts allows separation of pressure peaks due to wave reflection from shock tube walls 
from the pressure history. These pressure spikes can then, for example, be removed 
from structural frequency excitation analyses. Therefore, the information obtained 
via numerical simulation can be used for design from dynamic similitude conditions, 
and to augment data obtained in the test facilities. 

The simulation of the blast-wave problem has been studied in varying degrees 
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of physical complexity, from self-similar Euler to two-dimensional viscous flows. The 
simulation of the blast- wave problem using the Euler equations was studied by Kutler, 
Sakell, and Aiello [2] in 1975, where the intersection of a planar shock with a wedge 
in supersonic flight was modelled using the MacCormack scheme [3]. The self-similar 
nature of these types of flow with respect to time was used to obtain two- and later 
three-dimensional cone solutions [4]. In addition, Kutler and Shankar [5] used a shock- 
fitting procedure for the regular diffraction of a planar shock by a wedge, which is 
also self-similar. Although good results can be obtained using discontinuity fitting 
methods, coding complexities have generally hindered their application for complex 
situations. 

A general solution of the truly time-dependent inviscid interaction problem was 
modelled by Champney, Chaussee, and Kutler [6] in 1982. Shock diffraction over sev- 
eral simple two- and three-dimensional geometries were presented using the 
MacCormack and Beam- Warming schemes [7]. Mark and Kutler [8] performed a 
two-dimensional simulation of a shock passing over a simplified profile of a truck. 
However, inaccuracies due to discontinuity smearing and oscillations led to the de- 
velopment and application of explicit and implicit high-resolution schemes in the 
mid- 1980 ’s. Several two-dimensional problems using these high- resolution methods 
were presented by Yee [9] and Hisley and Molvik [10]. Recently, Lohner [11] has ob- 
tained solutions for complex geometries via an adaptive unstructured approach. The 
geometric flexibility of an unstructured grid method offers great potential, however 
the use of stretched tetrahedral grids required for efficient computation of viscous 
flows is a current topic of research. 

The solution of the Navier-Stokes equations in the high Reynolds number regime 
requires the use of fine grids to resolve thin viscous layers. The concomitant stiff- 
ness arising from the difference in length scales suggests the use of implicit schemes. 
Bennett, Abbett, and Wolf [12] applied a Beam- Warming scheme [7] with the Baldwin- 
Lomax [13] turbulence model to the problems of a developing boundary-layer behind 
a moving shock and shock diffraction over a cylinder. Molvik [14] used implicit 
high- resolution methods to obtain two-dimensional solutions for the unsteady devel- 
opment of a boundary-layer, the cylinder diffraction problem, and intersection of a 
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planar shock with a missile in supersonic flight. 

The purpose of this effort is to address the general three-dimensional, viscous 
blast- wave problem. The techniques developed here utilize total variation diminishing 
(TVD) upwind and upwind-biased schemes to resolve the discontinuous flow features 
without the oscillations prevalent in the more conventional central difference methods. 
Wave speeds are resolved adequately at large Courant numbers through the use of 
time conservative differencing and Newton subiterations. 

1.1.2 Cavity Flow Effort 

The Stratospheric Observatory For Infrared Astronomy (SOFIA) will be a three meter 
class Cassegrain telescope which utilizes a Boeing 747SP as an observation platform. 
An artist’s concept of the observatory, which is a follow-on to the Kuiper Airborne 
Observatory, is shown in Fig. 2. This airborne system, currently being studied by the 



Figure 2: Artist’s concept of the SOFIA configuration 


United States’ NASA and Germany’s DARA, offers capabilities which augment land 
and space-based options in several ways. First, the mission flexibility of a long-range 
mobile platform lends astronomers freedom to investigate transient astronomical phe- 
nomena on a global basis. Second, atmospheric attenuation of some wavelengths 
of interest provide motivation for a platform which operates above the tropopause. 
Third, the cost of maintaining and upgrading observation technologies is lower than 
would be incurred with an orbiting configuration. 
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Nevertheless, the use of an aircraft-based observatory presents some challenges. 
The limited bandwidth of solid materials in the infrared frequency ranges of inter- 
est preclude their use as windows. The temperature of the window material would 
also contribute to background radiation levels. Therefore, the telescope cavity must 
remain open to the freestream. Empirical evidence has shown that violent shear 
layer oscillations with concomitantly dangerous levels of acoustic loading occur for 
untreated, or rectangularly shaped, open cavity configurations [15, 16]. Hence, there 
is a need to develop cavity flow control treatments to suppress the flow unsteadiness, 
both to reduce the risk of injury to the crew and to obtain high-quality seeing. To- 
wards these objectives, both experimental and computational fluid dynamics (EFD 
and CFD) analyses will be used in the design cycle. The purpose of this work is 

to develop and apply numerical tools for use in the design of the next generation 
airborne observatory. 

The driven cavity problem has been a subject of much research, both experimental 
[15-29] and numerical [30-40], owing to its wide practical applicability. The buffeting 
and sound production of bomb bays, slotted wind-tunnel walls, transition-delaying 
airfoil cavities, and deflected control surfaces are examples of the range of problems. 
The effort here is focused upon the transonic regime, and previous efforts in these 
flow speeds are reviewed by Komerath, Ajuha, and Chambers [41] and Rockwell and 
Naudascher [42]. Some of the research which is pertinent to the transonic aero- window 
problem is highlighted below. 

In 1955, Karamcheti [17] studied subsonic and low supersonic flow over rectangular 
cavities, in which the inverse relationship between cavity length, L, and dominant res- 
onant frequency, as well the acoustic intensity and radiation directivity was reported. 
In the same year, Roshko [18] documented skin friction and pressure distribution along 
the cavity walls. The three-dimensional low subsonic study of Maull and East [20] 
showed that spanwise cells can modestly perturb C p distributions from the idealized 
infinitely wide cavity. Rossiter [15], in 1964, related cavity resonance to the edge-tone 
phenomenon, and deduced a model applicable to the transonic regime of concern here. 

In this study, Rossiter also demonstrated that a cavity leading edge spoiler drasti- 
cally reduced acoustic levels. The work of Heller and Bliss [24] demonstrated that 
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the cavity feedback mechanism can also be suppressed by geometry modifications at 
the shear layer impingement region. These wind tunnel tests have provided valuable 
insight into the governing flow mechanisms, but dynamic similitude was typically not 
achieved. The dangerous behavior of cavity flows generally precludes the use of flight 
tests to verify scaling laws. Numerical simulation offers the potential for safely scaling 
sub-scale tests to flight conditions. 

Toward this objective, numerical efforts of the type which can model the viscous 
flowfield about a geometrically complicated structure was begun in 1979 by Hankey 
and Shang [30]. Using MacCormack’s scheme, the dominant resonant mode of a rect- 
angular cavity was accurately predicted. Om [37] used the same scheme to compute 
flow about quieted two-dimensional cavities. In 1987, Suhs [34] used a block im- 
plicit scheme to obtain the viscous flow about a parallelepiped cutout. The overset 
mesh method which was the precursor to that used herein was utilized. Dougherty 
et al. [39] have recently completed a detailed study of two-dimensional cavities using 
a high resolution scheme. They computed distinct spectral peaks which have been 
observed experimentally. 

The present effort builds on past numerical studies by validating the ability to 
predict free shear layers, and both untreated and treated two- and three-dimensional 
cavity configurations with an efficient scheme in an overset mesh framework. Compar- 
ison of computed flowfields to experimental and analytical results allows assessment 
of cavity load prediction capabilities. Prediction of the acoustic intensity levels and 
frequencies are of primary interest for safety reasons. However, estimation of optical 
distortion is required to determine mission effectiveness. 

1.1.3 Aero-Optics Work 

The study of the effect of a fluid field upon an optical field, dubbed aero-optics, has 
been extensive over the past four decades. Applications include imaging of re-entry 
vehicles or, as in this study, of astronomical bodies through the atmosphere. Many 
experimental and theoretical approaches to the optical distortion problem have been 
investigated, those of which are pertinent to this transonic aero- window problem are 
summarized here. The experimental efforts can be grouped into two categories: direct 
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measurement methods and techniques based on aerodynamically inferred quantities. 
Results obtained via the latter method are more prevalent because of practical diffi- 
culties in direct measurement techniques [43], In fact, only aerodynamically inferred 
distortion levels will used for validation in the present work. 

Although early experimental and theoretical efforts assumed incoherent statistical 
turbulence [44, 45, 46], recent studies have begun to examine the effect of shear layer 
structures on electromagnetic field distortion. Using a passive scalar field from a 
direct numerical simulation, Truman and Lee [47] found an optimum viewing angle 
normal to the hairpin vortices in the homogeneous sheared fluid region. They also 
found analysis via non-refracting geometric optics to be equivalent to the parabolized 
Helmholtz representation of light. Although this class of studies provides excellent 
insight into the effects of small-scale structure on the electromagnetic field, it is clear 

that the expense of such methods precludes their near-term use for the problems 
under consideration here. 

The study of large scale structures in shear layers has been an active topic of 
research since they were observed by Brown and Roshko in 1974 [48], Only recently 
has the effect of these structures on the optical field been studied. In 1990, Chew and 
Christiansen [49, 50] experimentally observed the effect of shear layer structures on 
beam propagation. Tsai and Christiansen [51] used an Euler simulation to determine 
the optical characteristics of a perturbed free shear layer. The use of a growing 
sinusoidal phase plate to represent the effect of vortical structures on an optical field 
was hypothesized. Wissler and Roshko [52] recently performed an experimental study 
of the motion of a thin light beam caused by passage through a shear layer. They 
postulated that spanwise steering asymptotes to a higher level than the streamwise 
component. 

The numerical modelling of the optical effect of a cavity-spanning shear layer was 
presented by Cassady et al. [53] in 1987. They found their two-dimensional solution 
to result in poor prediction of optical distortion. Farris and Clark [54, 55] used time- 
mean quantities and empirical evidence to ascertain the fluctuating density levels 
required for optical analysis. 
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The present effort attempts to determine what portion of the optical path distor- 
tion can be resolved using cell sizes required to obtain an accurate flowfield solution. 
Computed optical distortion levels are compared to flight and wind tunnel measure- 
ments for two- and three-dimensional quieted cavities, respectively. The following 
chapters address the methods used to predict the unsteady flows, the modelling of 
turbulence, and the optical distortion. 



Chapter 2 

Numerical Method 


2.1 The Governing Equations 

The Navier-Stokes equations may be expressed in integral conservation law form, 
coupled with the continuity and energy equations as 

Uvl Qdv ) + vi pds=0 (d 

where body forces have been neglected and the cell volumes are time invariant. Here 
V is the volume of an arbitrary fluid packet, F = i?jvsi + ^/vsj + G/vsk is the flux 
tensor of second order, and ds is an outward directed normal of a differential surface 
area. The vectors may be written in Cartesian coordinates as 


Q = \p, pu, pv, pw, ef 


pu 

PU 2 + p + T xx 

EnS = | puv + T xy 

PUW + T xz 

(e + p)u + T xx u + T X yV + T xz w + q x 


9 
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pv 

PVU + Ty X 

FnS = pv 2 +P+ Tyy 

PVW + Ty Z 

(e + p)v + Ty Z U + TyyV + TyfW A 

pw 

pwu + T zx 

GtfS = pwv + T zy 

pw 2 + p + T zz 

(e + p)w + T ZX U + T Z yV + T ZZ W + q z 
where each flux can be partitioned into inviscid and viscous portions. The density, 
pressure, and velocity components are respectively given by p,p,u,v, and w. The 
viscous stresses are composed of the terms: 


T xx 

T VV 

T xx 


. du ( du dv 

— 2P"X A I — — + — h 

p dx l dx dy 



-2p— - A 

dy 


du dv dw 
dx + dy+ dz 



{ du dv dw\ 
+ dy + dz ) 


f du dv\ 

Txy = Tyx ~~^{d^ + di) 

( du dw\ 

Txz = Tzx = ~ ti [d^ + d^J 

( dv dw\ 

r y , = T„ + 


The total energy per unit volume, e, is related to the internal energy per unit mass, 
e, by e = pe + pq 2 /2. The perfect gas equation of state, p = pRT, completes the 
system. In addition, for thermally and calorically perfect gases, the internal energy 
per unit mass and the enthalpy per unit mass can be expressed solely as functions of 
temperature: 


de = c v dT, dh = CpdT 
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Finally, for a calorically perfect gas the specific heats are constant, leaving e = c v T + 
const., and h — c p T + const., where the additive constants may be set to zero. The 
ratio of specific heats and specific gas constant are 



= - for air, R = c p -c v 


and the thermodynamic variables are related using 


h T= (i±P}' e = -^. + |(„2 + 


where h j* is the total enthalpy per unit mass. 

Fourier’s law for heat transfer by conduction is assumed; hence, the heat transfer 
can be expressed as 


q — kVT — - ( q x \ + q y j + 9 z k) 
fde.de. de\ 

~ K W + d^ J + dI k ) 


where k - k/c v - m/Pr. The Prandtl number for air, which is a function only of 

the gas, relates the diffusion of momentum to the diffusion of heat, and is fixed at 
Pr = 0.72. 


The relationship between the first (/z), second (A), and bulk (C) viscosity coeffi- 
cients is C = in + A. The bulk viscosity coefficient is set to zero in accordance with 
Stokes’ hypothesis, resulting in A = -§^ This hypothesis is invoked here based on 
the assumption that the relative effects of the shearing stress is much larger than those 
caused by the dilational stress effects, not on the theory for monatomic gases [56], 
By using the kinetic theory of gases, the physical phenomena of thermal conductivity 
and viscosity can be expressed in terms of the thermodynamic states. Viscosity is 
related to the thermodynamic state using Sutherland’s formula: 


where C\ and C^ are specific to the 


C{T\ 

n = 

t + c 2 

gas in question. 
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2.2 Turbulence Model 

Current limitations in computing power relegate most engineering computations to 
the use of grids which are too coarse to resolve all of the pertinent scales of motion. 
Reynolds averaging of the Navier- Stokes equations decomposes the flow into slowly 
and rapidly varying components [57]. The slowly varying component is resolved from 
the spatial and time integration step sizes, while the rapidly fluctuating component 
is modelled. Although the blast-wave computations did not model turbulence, the 
turbulence model used in the cavity flow problems is outlined here for completeness. 

The effective or eddy viscosity due to additional turbulent mixing can be related 
to mean stresses using the Boussinesq approximation. The total effective viscosity is 
then given by ptotai = Pmoiecuiar + ^turbulent = P + Pt- The Reynolds stress resulting 
from the Boussinesq assumption is 



where ( ) denotes a time mean. The eddy viscosity is given by p t oc piv where the 
length and velocity scales are given by i and v. Alternatively, the expression for the 
eddy viscosity given by Prandtl is Pt oc pf 2 ]^, where the magnitude of vorticity is 



The local density is specified by Morkovin’s hypothesis, which states that compress- 
ibility does not affect the scales of the turbulent motion. 

The algebraic turbulence model of Baldwin and Lomax [13], as modified and 
implemented by Buning [58], is described below. This description is included to clearly 
show how the modified model constant used in the computed cavity shear layers is 
determined. The description assumes flow in the ( x,y ) plane with the freestream 
aligned with the x coordinate. 

Treatment of Wall Bounded Flows 

In Prandtl’s mixing length hypothesis the turbulent eddy size is limited by the prox- 
imity to the wall, giving t — ky, where the von Karman constant k = 0.4. The 
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addition of Van Driest wall damping results in 


£ = ky (l — e y+ ^ + j 

where A + = 26, y + = s ' P ^“’ y , and the wall shear stress is t w = • 

The Baldwin-Lomax two-layer model uses the Prandtl-Van Driest model for the 
inner layer, and is given by 


Pt 


pt 2 M y < ycrosaover Prandtl-Van Driest 

pCcpH EwakepKleb V ^ V crossover Outer region 


where y is the distance from the wall, ycrossover is the location of the first intersection 
of inner and outer values of fi t , the Clauser constant is K = 0.0168, = 1.6, and 


F wa k e = min 


( 


ymaxF rr 


C wkUn 



where C w k = 1.0. The quantities F max and its location y max are found from 


F (») = M (l - 


where the exponential term is dropped in wakes. The search for the F max term ends 
when F(y) drops below a specified percent of the first peak away from the wall, or at 
overset mesh boundaries. The Klebanoff intermittency function is specified according 
to 

Fxuiv) = 

where Cxub = 0.3. The total velocity difference is given by 



Udif = \/(« 2 + v 2 ) - \/(u 2 + v 2 ) f 


and the latter term is taken as zero at the wall. 


Treatment of Free Shear Flows 

The eddy viscosity in the shear layer was computed as outlined by Buning [58]. 
Development of the free shear layer model begins by using F(y) = y\u\, as suggested 
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by Baldwin and Lomax [13] for wake regions. This results in 


F wake — C w lc 


ymaz^dif 


= c.tfrr*-) m— 

VMma*y 

oc £ 2 |u>| 

where specification of C w k is discussed below and the velocity difference is modified 
to be half the total velocity difference between the streams in the specified shear layer 
region 

“*/ = v/(« 2 + " 2 )„„ - >/(«’ + 

Finally, the Klebanoff intermittency function is modified to 

( CjOet\y - j/|u,[ r 


F K ,eb(y) = 


1+5.5 


V 


Vw 


! ) 


where the shear layer width is given by y w = u dif / |u>| roax . The free shear layer model 
is now given by 

n t = pK CcpC W k r~r^— (2) 


M 


after dropping the intermittency function for the analysis below. 

The magnitude of the eddy viscosity in the free shear layer model can be altered 
by specification of C w k • The remainder of this section shows how C w k is chosen based 
on empirical and analytic information. 

Gortler’s shear layer solution is given by 

ui + u 2 f, , u 2 - ui ,1 2 f x _„j , ay 

u = — — [1 + /(,)j . erm = ^ /„ « ’ *», 1 - T 

where U\ and u 2 are the velocities of the slow and fast streams and 77 is the similarity 
coordinate (See Fig. 19). The spreading parameter a is inversely related to the 
spreading rate, db/dx, where b is a measure of the shear layer width. The value of 
the spreading parameter when the velocity of one of the streams is zero is <tq. 
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Gortler s solution can be used to determine the maximum vorticity magnitude in 
a free shear layer as follows: 


U) I = 

du 


du drj 

' 1 

dy 


dydy 


a j 

Aue-" 


M 


Xy/ft 

° A 
Am 


, A U = («2 - Mi) 


XyfH 


Now, using Prandtl’s mixing length assumption and scaling laws for jet bound- 
aries, eddy viscosity can also be expressed as 


Ht °c pl 2 M 

2 Am 

K pl T 

= KopbAu 


( 3 ) 


where K 0 = 


46^? 


Setting Eqs. 2 and 3 equal results in 


C wk = ( a 0 KC cpy/t)- 1 

and only a 0 remains to be specified. Estimates of a 0 from empirical evidence is quite 
variable, ranging from 9.0 to 13.5 primarily dependent upon whether the upstream 
boundary layer is turbulent or laminar [48, 59, 60]. For this series of cavity flow 
efforts a 0 was set to 11.0, which appears to be the result from the highest quality 
experiments, resulting in a value of C wk = 1.91. Previous numerical investigations 
appear to indicate that capture of resonance is not strongly dependent upon the 
turbulence model in the cavity [30, 34, 36]. 


2.3 Transformation to Curvilinear Coordinates 

In order to adequately resolve the solid boundary/fluid interaction, it is common to 
transform the governing equations into curvilinear coordinates which can be body- 
conformal. Specifically, the body is constrained to lie at a constant f , 77, or £ level. 



CHAPTER 2 . NUMERICAL METHOD 


16 


For a stationary grid, this transformation can be expressed as 

r = t, Z = Z(x,y,z), 1 = T](x,y,z), C = C {x,y,z) 

Application of the chain rule of differentiation yields 

d d 

di~ u d( ’’'a^ + tl 9c 

with similar expressions for the partials with respect to y and z. The inverse trans- 
formation gives d d d d 

F( = H di + V( ai + Z( e~z 

Again, expressions can be found for the 77 and C partials in a like manner. Represented 
in matrix form: 


d_ * 
dx 


'Cr 

lx 

Cx ' 



d_ 

dy 

= 


Vy 

Cv 


d_ 
d 77 

f 

8>l® 

S 

.6 

lz 

G . 

✓ 

d 

dC J 


and for the inverse transformation, 


d I 


- -| 


9 1 

dZ 


y« H 


dx 

d 




d_ 

d 77 


■Xr) Vi) ■Zf 1 


dy 

d 




d_ 

. dC . 

\ 

. x < y< z < . 

✓ 

. dz - 


T-» 

Combining the use of T = (T -1 ) -1 and finite volume metrics, such as those described 
by Vinokur [ 61 ], leads to a scheme which is freestream- preserving because of the 
telescoping property. Hence, if the surface normals to a constant £,77, or C plane are 
defined respectively as 

s i+ l = Sx.i+1 1 + S y,i+fj + S *,<+i k 

= i(r 7 - r 4 ) X (r 8 - r 3 ) 
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s j+j — 

= ^( r 7 - r 2 ) X (r 3 - r 6 ) 

S *+l = 5 x,*r+lj + 5 J/,/t+lj + S z,k+$^ 
= ^( r 6 - r 8 ) x (r 5 - r 7 ) 
where the index convention is shown in Fig. 3. 



Figure 3: Hexahedral cell and stencil 
The metrics can then be formed as 

t* = JivnH-v&n) = ^s xi+ 1 

Cv = J(x XZn -***<) = y s y,i+i 

6 = J{x^ - x<yj = is Z)l+ 1 

rj z = J{y<n ~ y { *<) = ^s xJ+ i 

*lv = ~ HH) = ^ s yj+i 

Vz = J{x c yt-x ( y c ) = ys, J+ ! 

Ci = J(yt z i) ~ yti z 0 — y5 I>jt+ i 
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Cy = J{ X T] Z t ~ x i z r)) ~ y 5 y,fc+l 

G = ~ X T ) Vi) = y S *,fc+l 

These metrics represent the projections of the cell face normal into (x,y, z) space. 
The faces of the hexahedron exactly enclose the discrete control volume, i.e., no gaps 
axe permitted at the edges. 

Finally, the Jacobian of the coordinate transformation is equivalent to the inverse 
of the volume, as related by 

1 _ d(x,y,z) 

J d(t,r),0 

= x ((yr, z c - y< z v) - x v(yt z < - y< z t) + x < - y^n) 

= V = -(s i _i+s i _i+s /t _i).(r 7 -r 1 ) 

Utilizing these metrics in the application of the chain rule to Eq. (1) and subse- 
quent simplification yields 

q' t + E'f. + f; + G" c = 0 

where 

Q' = QV 

E'ns = (Ens€x + Fp/s^y + Gfifstz) V 
= [EffS s x + Ffl/sSy + 

F'nS = (FnSVx + Firstly + GnS^z) V 

= + FffsSy + Cxjvss*].? 

G'nS = (EnsCx + FffsCy + GfisC) V 
= + FffS s y + GffsS z ]k 


Separating the inviscid and viscous portions of the flux vectors, then in the f direction 
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E' ns — E’ + E[„ where 

pU 

puU + £ x p 
E' = V pvU + (, y p 
pwU + £ z p 
L (e + p)U . 

Here the contravariant velocity component in £ is U = u^+v^ + tv^ without metric 
normalization. The viscous flux can be represented as 


0 

1~xx£x d" T yx£y “I" Tzx£z 

T xy£x + Tyy£y + T 2 y£ z 

1~xz£x 4* Tyz£y + T zz £z 

(ue' 2 + ve f z + we 4 ) + ( q x £ x + q y £ y + q z £ z ) 
where the viscous stress terms are evaluated by again invoking the chain rule, and 
the flux in the 77 and C directions are found similarly. The results presented herein 
axe implemented using either the thin-layer or the full viscous term treatment. 

The widespread use of the thin-layer approximation, first implemented by 
Steger [62], can be justified from either physical or algorithmic arguments. Physi- 
cally, the neglect of all diffusion processes parallel to the body is similar to that used 
in boundary-layer theory, albeit not as restrictive. Hence, when the viscous effects 
axe confined to thin regions along a constant £, r?, or C plane, this assumption is valid. 
Regarding the algorithmic argument, the banded matrix structure used in multidi- 
mensional algorithms which sequentially solve a set of unidirectional problems can 
include only these thin-layer terms implicitly. This thin-layer flux in the 77 direction, 
assumed to be the body normal coordinate, is expressed as: 

0 

m iu 7 , + m 4 v v + m^Wr, 
p> — _y m^Ujj + m 2 v n + m 6 w v 

m^Ur, + m 6 t;, + 

+ m 2 vv r , + m^ww^ + m 4 (uv v + vu v ) + m 5 (u^ + w Ur) ) 
+m 6 (vw r , + u>v n ) + Ke^l + ^ + *7?) 
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where the (~) denotes an arithmetic mean value and 


m\ 

= p\ 

+ V 2 v + »?z) 

, m 4 = 

P 

g r lx r ly 

m 2 

= p i 

[vl + \nl + vi) 

, = 

P 

m 3 

= p i 

(*jj + i 6+ \vi) 

, m 6 = 

P 

z w 


These viscous flux terms may be found for the remaining spatial coordinates as well. 
The results presented here are implemented using either the thin-layer or the full 
viscous term treatment, as required by the the flow physics. 


N ondimensionalizat ion 

The governing equations may be nondimensionalized by the choice of a length scale, 
denoted by Z, and reference values of p, u, and p such as 

Pref = Pool Uref = yPoo/ Poo 5 Pref — Poo 

The nondimensionalized variables follow: 

P = p/pref, P = PlPref, C = e/(p re fU 2 ref ) 

U = u/iiref, V = v/u re f, W = w/u Te f 

t = tU re f/L, T = p/p, P — pf Pref 

The Reynolds number resulting from this procedure is Re = p re jLu r ef I Prefi where 
the (“) denotes a dimensional quantity, and ( )«, denotes the freestream conditions. 


2.4 Upwind Schemes 


The numerical scheme will be described using first-order terms, following which the 
higher-order extensions will be outlined. The scheme expressed for a cell which has 
a mean flux value on each of the six sides is 


7rJ QiV 

or Jv 
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+ 

+ 



- 

- GU k _^)dT)di = 0 


In fully discrete form, after dropping the primes for convenience, the governing equa- 
tions can be written as 


+ 

+ 



= 0 


where n denotes the time level in this implicit representation, and A£, A 77 , and A( 
are set to unity for convenience. 

These flux terms may be evaluated using a technique which may be broadly classed 
as either central or upwind. The latter technique is chosen for this study for the de- 
sirable numerical properties, such as diagonal dominance of the flux Jacobian, and for 
the physical dependence on zones of influence which are inherent in upwind schemes. 

Upwind schemes bias the derivative evaluations required to determine the flux 
across fluid cells according to the sign of the characteristic speeds. In this manner 
these methods bring the physics of the hyperbolic system, the unsteady Euler equa- 
tions, into the numerical solution process. To facilitate the implementation of these 
upwind schemes, the eigensystem is determined. The similarity transformation which 
diagonalizes the unsteady, inviscid, gas-dynamic equations, shown by Warming, et 
al. [63], is outlined as follows 




dQ 


= TAT" 


where the rows of T 1 are the eigenvectors and 


A = A + + A - = diag[0,U,U,U + c,U -c] ||s|| 

using normalized contravaxiant velocity components. The eigenvalues can be split 
according to, among other splittings, their signs: 

A±|>| 


2 
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where A is an element of A. 

Two upwind schemes are implemented here to compare the results which may 
be obtained with either of the techniques. The initial portion of this discussion will 
be presented unidimensionally for simplicity; the multidimensional extension will be 
outlined towards the end of this section. 

2.4.1 Flux Vector Splitting 

The shock-capturing scheme developed by Steger and Warming [64] revisited the 
classical characteristic procedures. They found that the Euler equations possessed the 
property of homogeneity of degree one for the equation of state used here, meaning 
E(aQ) = aE(Q). For a vector with this property E = AQ, where A is the flux 
Jacobian given by dE/dQ. Consequently, the flux vector can be split into two parts, 
each physically corresponding to the right and left moving waves. This technique 
resulted in the flux being represented as a combination of the subspaces associated 
with the positive and negative eigenvalues, expressed as 

E = 7 , (A + +A-)T- 1 Q = (A + + A~)Q 
= E+ + E~ 

where T and T -1 are the right and left eigenvectors of the flux Jacobian matrix A , 
respectively. The flux across a cell face can be determined by 



Because the Jacobian at i + dependent on two states, this solution method now 
diverges from the original Steger- Warming flux vector splitting. The treatment of 
this Jacobian is shown in a following section. Linearization in time can be performed 
in one dimension as follows, extension to the multidimensions is straightforward and 
is omitted for brevity. At the n + 1 time level, 


£,"+} = (A+)’»Qr l + (A-^egi 1 

= M + )r4«3r+(^')rH*o? + i+ B r + j 
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where the implicit change in the dependent variables is given by SQ n = Q n+l — Q n 
Note that the Jacobian matrices are frozen at time level n. The remaining flux, E n+ } , 
may be obtained similarly. ' 2 

To assess the effect on stability of this type of linearization, a procedure developed 

by Barth [65] is applied to this method. Using the semidiscrete form dQi/dt = — i? t , 
then 

Using frozen Jacobian matrices, the method can be linearized as follows: 

[s + (§)’)«—*■ «> 

where for the first-order subset the Jacobian is a block tridiagonal matrix. The blocks 
along the i th row are 

— = —(-A* \ 
dQi-i Ax V *-j/ 

— = — N + -a- \ 

dQt Air l '+i -i) 

— = —(a- ) 

dQi+i Ax \ ,+ 2 / 

This scheme is inherently conservative in space because of the telescoping property of 
the finite- volume formulation; analysis of this scheme reveals that it is also conserva- 
tive in time, possibly allowing the use of large Courant numbers [65]. A demonstration 
of this analysis proceeds by writing the scheme as 

[^+^ n ] SQ = -R n = -A n Q n 

hence 

Q n+1 + AtA n Q n+1 = Q n 

In order for the scheme to be conservative in time over a periodic domain, the global 
average of a solution must remain constant for all time, i.e., £f-i ^ t Q n = 

£i=i Q? +1 - Hence, when Eq. (5) is summed across the domain, the result is that 
the columns of A n must sum to zero. For example, summation across a three-point 
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domain yields 


[Qi + Q 2 + Q 3 ] n+1 +At[7, /, /] 

' b \ c; . ' 
a; b \ c; 

n 

Qi ' 
Q 2 


. a ] b ; 


.<? 3 . 


[Qi + Q 2 + C?3] n 


where the elements of the middle column of A n are C[ = A y 2 , B 2 = A^ 2 ~ A 3 / 2 , 
A” 3 = —^ 5 / 2 ’ which sum t° zero. 


2.4.2 Flux Difference Splitting 

Flux difference splitting methods are based on the Riemann problem, solved exactly 
by Godunov in 1959 [66]. The Riemann problem is composed of m + 1 piecewise 
constant states separated by m wave families. The waves include shocks, contact 
surfaces, and rarefaction fans. For each of the Riemann problem cells, the transition 
of the dependent variables is a function of a parameter family. The solution can be 
found once these transition states are known. Approximate Riemann solvers simplify 
the numerics of the problem by eliminating the iterative process required to find the 
intermediate states. 

t 



Figure 4: Riemann problem schematic 


Figure 4 shows a schematic of the Riemann problem with the piecewise constant 
states separated by the appropriate wave families. The flux through the cell face is 

1 


E«i = 


2Ax 

1 


B, + £, + ,-(A£^-A£- i )] 
2Ax (£i + Ei+1 _ A|£|i+ i ) 
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where |i?| — \A\Q = (A + — A )Q = [T(A + — A )T l \Q. The flux differences associ- 
ated with the + and — traveling waves are 


AE+, = (rA + T-') j+ j(Q j+1 - <?,) 

AE-J = (rA-r-‘) i+J (Q i+ , - (?,) 

Again utilizing the semidiscrete form dQifdt = — i2,, then 

= 55; [ (£,+i " - (A|E| -+i + A i £ i<-t>] 

This method can be linearized using a procedure similar to that described previously 
in Eq. (4). The blocks of the i th row are expressed as 


dR, 

dQi-i 

dR, 

dQi 

dRi 


1 

2Ax 

1 

2Ax 

1 


—Ai-i + 


dA\E\i- 

dQ 




aAiEij.i 


A {+ 1 — 


dQi ' dQ 
dA\E\ i+ 






dQ i+ 1 2Ax V"” ri dQ i+1 
Substitution of the flux difference splitting expression yields 


5A|£|,_ 


dQ< 


t+1 


dQi+i 


- Mlt+i + 


[Wi+jWi+l - «•)] 

dQi+i 


(Qi+i - Qi) 


These true Jacobians are expensive to compute, and the simplification to approximate 
Jacobians is made as 


gAgUi 


dQi 


i+l 


MI.+4 


Utilizing these approximate Jacobians, the linearization proceeds as 


P ' n + 1 — 1 

2Ax 
1 


2Ax 


(£r +1 + ET+Y - A|E|" + + i) 

,T 2 

[W + + (AJ., - W +1 )«? i+ i] + E" + 1 


where E n+1 = E n + A n SQ and Q n+1 = Q n + ( dQ/dt) n At . The fluxes through the 
remaining faces are determined similarly. This scheme can also be shown to obey the 
criterion for conservation in time. 
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2.4.3 Roe Averaging 

In order to determine the Jacobian at the cell face i + i, some function, A = 
MQl, Qr), must be assumed, where the subscripts indicate left and right states. The 
location of this flux evaluation is one of the differences between finite-difference and 
finite- volume schemes. The evaluation used here is attributable to Roe [67], which 
provides an approximate solution to the Riemann problem. This Jacobian is cre- 
ated through the use of a parameter vector composed of a geometric-like mean of the 
states. The more obvious arithmetic Jacobian forms, such as A = \(Ai + Ar) or A = 
A(^(Ql 4- Qr)), are not conservative forms. Conservative Jacobian forms satisfy 
A{Ql, Qr)(Ql ~ Qr) = Ei — Er. Stated explicitly, the Roe averaging operation is 

P — \JPlPr 
- _ m «£\/pI + UjRy/PR 

yfpi + y/PR 

_ PL^it + pjUii + Uj R ) + PRUjR 

Pl + 2(> + Pr 

^ _ PLjJjj^L + p((^r)i + ( hr);? ) 4- PR^ha^R 

pL + 2p + Pr 

where a (“) denotes a Roe averaged quantity, and the latter forms are presented as 
inexpensive alternative expressions. Substitution of (i) for L and (i + 1) for R allows 
the evaluation of the Jacobian at the intermediary cell face. For the flux vector 
splitting case described earlier, MacCormack [57] has found this average helps to 
alleviate excessive numerical dissipation in regions dominated by viscous effects. Roe 
averaged values axe utilized throughout the development presented here. 

2.4.4 Higher-Order Extensions 

Spatially first-order methods frequently provide inadequate resolution of the flowfield. 
However, the methods discussed above can be extended to higher-order spatial accu- 
racy by modification of the right hand side. In order to assist in the preservation of 
well-behaved solutions near the discontinuities admitted by the strong conservation 
law form of the Euler equations, a total variation diminishing technique is imple- 
mented. If the total variation of a solution is defined as TV(u) = |u,+i — u,|, 
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then a solution which follows TV(u n+1 ) < TV{u n ) is TVD. The TVD constraint can 
be shown to result in diagonal dominance, allowing the use of relaxation schemes. 
In this manner the scheme may be extended to higher space accuracy throughout 
the smoothly varying regions of the field, reducing the accuracy in localities of high- 
gradient and extrema in order to obtain sharp and oscillation-free resolution. These 
methods are rigorously applicable only to scalar nonlinear equations or a system of 
linear equations in one spatial dimension. Application of these schemes to multi- 
dimensional systems of nonlinear equations are generally not TVD. Moreover, it is 
not clear that the higher-order accuracy of the unidimensional problem is retained in 

multidimensional cases. However, the results which can be obtained demonstrate the 
usefulness of the technique. 

Of the several methods which fall into the TVD domain [9], the technique im- 
plemented here is one attributable to Chakravarthy and Osher [68], the development 
of which follows for completeness. In this formulation, the higher-order flux can be 
expressed as a sum of a first-order flux, denoted E i+ i, and a flux correction term. The 
flux correction terms are determined by first computing the flux differences across the 
m wave families mentioned previously. Subsequent limiting of these flux differences 

and summation across the wave families results in the higher-order flux. This flux is 
expressed as 


E H 



+ 


(1-0) 

* 


(1 + 0) 
4 

E dE ‘H 

(1+*) 

4 

[fs&j 

+ (l-« 

4 

t • 

r m ~ J+ 1 

\Z iE <-i 
- * 


where (~) and ( = ) indicate a quantity that has been limited, j is the index denoting 
the wave family, and i is the index assigned to a cell center. Using the notation of V 
for the rows of the left eigenvector matrix, T~\ and for the columns of the right 
eigenvector matrix, T, then the measure of the change in the dependent variables is 


"i+i ~ l U\ ' (^h- 1 - Qi) 

The measure of the change in the flux is defined as 


a j = AV = (A i+ + \ j ~)oP 
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the eigenvalues being split as shown previously. The limited counterparts of these 
values are obtained as: 

= minmodfa-j./Jff'J 

d~ +h = minmod[<7- + i,/?(7 t - f ] 

= minmod 

K-i = minmod [a + _i,/?a + + i] 

This limiter returns the argument of smaller magnitude when the signs are equal, and 
returns zero when the arguments are of opposite sign. This procedure effectively adds 
dissipation locally in regions of high flux gradient and at inflection points. In this 
manner, monotonicity is preserved by preventing the creation of new extrema while 
preserving the global accuracy of the solution. While formal accuracy estimates are 
difficult to ascertain because of the nonlinear application of limiting to different wave 
families, numerical experiments have demonstrated that the global accuracy of the 
underlying scheme is preserved [69] . 

The compression parameter, /?, is restricted according to 1 < 0 < frf and tlie 
limiting operator is given as 

minmod(x, y) = sign(x) (max{0,min[|x|,ysign(x)]}) 

The compression parameter reduces the amount of dissipation added, the range being 
bounded by accuracy and TVD constraints. Finally, the limited flux difference values 

are expressed as 



= 

dE i+$ 

= ®!+} r .+s 

ML, 


dE { _ i 
1 2 

~j+ 

1 2 1 2 


This asymmetric limiter is designed to modify the fluxes only in the rapidly varying 
portions of the flow, where nonphysical oscillations are likely to occur. Since these 
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high-gradient regions are confined to thin regions, the dominant solution domain is 
differenced in accordance with the underlying scheme. Variances in the value of the 
compression parameter allow the fluxes to be limited for different gradient levels. This 
implies that use of 0 max will cause the limiting action to be taken only in the high- 
gradient regions, and lower values of 0 will result in limiting for commensurately lower 

flux gradients. The variety of schemes which can be obtained using this technique 
are shown in Table 1 . 


4> 

Unlimited Scheme 

Pmax 

2 nd order TE 

-1 

1 

3 

0 

1 

3 

1 

2 

1 

Fully upwind 

Nameless 

Fromm’s 

3 rd Order 

Low TE 2 nd Order 

Central 

2 

5 

2 

3 

4 

5 

oo 

K A *) 2 /~* 

g(A xff xxx 
n(A x) 2 f xxx 
0 

-i(A X) 2 f xxx 

~U Ax ) 2 f*** 


Table 1: Summary of schemes 


Here TE = (5 — <f>){Ax) 2 f xxx /4 defines the leading term of the truncation error 
for the unlimited form of the schemes. Local metrics have been used in the above 
method to maintain reasonable computational efficiency, a satisfactory approximation 
for grids which do not contain rapid variations. 


2.4.5 Viscous Terms 

The viscous terms are treated through central differencing about the cell faces. The 
explicit terms are conventionally differenced after chain-rule expansion, inclusive of 
the cross terms if these diffusion processes are significant for the problem at hand. The 
left hand side does not include these cross terms, and the resultant viscous Jacobian, 
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employing V = [p, u, v, w, e] T as the primitive variable vector, is 


M = 



0 

0 

0 

0 

0 

0 

m 2 2 

PS x&y j 3 

ps x sj 3 

0 

0 

m 23 

77733 

PSyS 2 j 3 

0 

0 

m 2 4 

77734 

77144 

0 

0 

77752 

77753 

77754 

77755 


m 2 2 = ri^ s l + + «z) » m 33 = K4 + + s]) 

m 4 4 =M 5 I + S y + 3 S z) > m 55 = K ( 5 x + S y + S z) 


m 52 = um 22 + wm 2 3 + wm 2A 
m 5 3 = tim 32 + u 77733 + 111^34 

77154 = UTT142 + 1777143 + Wm 4 4 

The viscous flux through a cell wall at j + \ is of the form Mj + i(Nj+iQj+i — NjQj) 
where 

0 0 0 ' 

0 0 0 

1 0 0 

0 1 0 

-v —w 1 

Now, using three-dimensional indices ( i,j,k ), expansion of the block structure gives 



N = 


dQ 

dV 


1 

P 


P 

—u 

—v 

—w 


0 

1 

0 

0 


I -e + k(u 2 4 - v 2 + w 2 ) -u - 
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+ 

+ 

+ 


+ B tj+k,k - B 7j- L,k 

+ c+ M 1 - c-' H 

^ V ^ij'Ic+L 

(yfe+i.fc " y M 'J+ | 6Qij+l ,k 

= AQiJ,k 


where only the thin-layer terms in rj are shown here. 


2.4.6 Factorization 


The extension of the techniques given above is accomplished through dimensional 
splitting. The method used here is that of Yanenko [70], where the factors are chosen 
in the £, 77 , and C directions. Expressing the three-dimensional equations in compact 
notation as 


<' + ti A+ f^ + s? c w=^ 


then the factorization procedure yields 


( i+ % A){,+ £, B)ii+ % c)6Q=AQ 


This system can be solved sequentially through the use of intermediary steps without 
loss of time accuracy. Although alternating direction implicit schemes of this type 
offer advantages of vectorization, the system is solved as a sequence of unidimensional 
problems, hence limiting the size of the time step due to stability restrictions [71]. 
The use of this technique here is justified by the requirement of adequate time history 
flow resolution, thus imposing an additional constraint on the maximum time step. 
Application of a line Gauss-Seidel method to the starting shock tunnel test case, 
discussed in Section 5.1.2, confirmed this hypothesis. This relaxation method offered 
a slightly increased stability range, but not enough to offset the additional expense 
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caused by short vector lengths for the shock tunnel problem. Additionally, since for 
the factored scheme the flux exchange occurs at the same time level, the technique is 
conservative, even when convergence at the subiteration level is not attained for each 
time step. 

The expense of both the upwind algorithms is relatively high: 86/zs per cell per 
iteration using a single processor on the Ames Research Center CCF Cray Y-MP / 832. 
These vectorized codes have computation rates of approximately 140 MFLOPS. In 
addition, the memory requirement is 40 words per cell. Decreased processor times 
may be achieved by many methods. For example, freezing the flux Jacobians for 
several subiterations will offer a processing time reduction of 15% per subiteration, 
albeit at the expense of memory. It is clear that the expense of these upwind methods 
is warranted only for problem classes in which the improved resolution is critical. 


2.4.7 Newton Iterative Technique 


Reduction of the linearization and factorization errors is achieved by a Newton iter- 
ative method of the type described by Rai and Chakravarthy [74] and Rogers and 
Kwak [75], albeit with the addition of allowance for a varying step size [14]. Assuming 
that the initial guess lies within the radius of convergence, the right hand side is con- 
verged to an arbitrary accuracy while holding time fixed. Since the right hand side 
includes the higher-order difference representations of the Navier-Stokes equations, 
linearization and factorization errors are eliminated at convergence. The method is 
discussed below where m is the Newton iteration index and n is the conventional 
index denoting time level. Discretizing Q t + E x = 0 gives 


l_(Qn+l,m+l _ Qn+l,m^ _ j. ^n+l.m+i _ Qn,m+\ _ (Qn+\,m _ 


« Q? +1 ’ m+1 - -}-(Q n+l ' m - Q n ) 
At 


where the solution is converged at time level n, hence Q n,m+1 = Q n . Defining bQ' — 

Qn+l,m+l _ Qn+l,m ? t h en 


I6Q' = ArQ" +1,m+1 — (Q n+1,m — Q n ) 

= -Ar£" +1 ’ m+ ' - (<2 n+1 ’ m - Q n ) 
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Linearization at iteration level m + 1 gives 

^jn+l,m+l _ ^£n+l,m _j_ 



n + ^m 


6Q' 


where the flux Jacobian has been frozen at iteration m. Substitution yields 
I6Q' = — Ar£’” +1,m - Ar-?-A n+l ' m 6Q' - (Q" +1 * m - Q n ) 


Rearranging results in 


Ax 


6Q' = — (Q" +1,m - Q n + Ar£'" +1 ’ rn ) 


which reverts to the standard noniterative form + ^j4"J 6Q' = —E£, when no 
subiterations are taken, as can be seen by substitution of n for n + 1, m. 

The temporally second-order accurate representation is found by extension of the 
above procedure. Using a three-point backward time stencil derived from a standard 
Taylor series approach, 


Qt = C 0 Q n+1 + C\Q n + C 2 Q 


n— 1 


where 


Co = 


C 2 = 


1 — o 


(1 - <t)At 2 -I- A-Ti 
-1 


, C, = 


(1 — <t)At 2 + Arj 


(1 - ct)At 2 + An 

and a — (1 + Ati/At 2 ) 2 . The elapsed time between the n — 1 and n time levels is 
given by Ati and between n and n + 1 is A r 2 . Rewriting at iteration level m + 1, 

C 0 (Q n+1 ’ m+1 - Q n+1 ' m ) = Q” +1 ’ m+1 - (C 0 Q n+1,TO + CiQ n + C 2 Q n ~ 1 ) 

Finally, 

\j + ——— 6Q' = — (Q"+ bm + £l Q n + C* Q n- 1) _ * £»+!,« 

L C/oaax j Co Co Co 

which reduces to 6Q’ = ^ (Q n -Q n ~ 1 )-E^ for the case of no iterations 

with fixed time step size. The formulation given above allows the use of a time step 
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size which is a function of time, but is fixed at each step for the entire domain. The use 
of a variable time step size allows the solution to progress using a constant Courant 
number, possibly preventing inadvertent divergences. Higher-order accuracy in time 
may be obtained by extension of the above technique, albeit with additional memory 
requirements. 

The assertion that this technique reduces the factorization and linearization errors 
is substantiated as follows. The right hand side of the method contains the discretized 
governing equations in their pure form, that is, without the numerical approximations 
utilized to attain rapid convergence. The left hand side allows the use of large time 
steps by relieving the Courant-Friedrichs-Lewy stability constraint. Deferring the 
question of uniqueness, if a set of dependent variables is found such that the right 
hand side is satisfied, then this field is a solution to the discretized equations regardless 
of the approximations made to arrive at that set. 


2.5 Diagonal Scheme 

Resolution of the transient wave-field of the blast- wave/target interaction problem 
class benefits from the use of upwind methods. In contrast, for the SOFIA effort a 
combination of the low transonic regime, the complex geometry, and the large prob- 
lem size required an efficient integration scheme. The algorithms used for the SOFIA 
effort, coded by Buning and Chan [58], are implemented within the Chimera over- 
set grid framework [72]. The solutions were obtained using a diagonal scheme [73], 
using spatially varying time steps for steady state computations, and fixed step size 
for unsteady flow simulations. The code utilizes the conventional dependent variable 
vector, Q = [p,pu,pv,pw,e] T , and pseudo-finite- volume metrics. Euler implicit time 
marching and second-order central spatial differencing were used for the computa- 
tions presented here. Computations were performed on the Numerical Aerodynamic 
Simulator (NAS) Cray Y-MP / 832 using SSD, at an expense of 14 ps per point per 
iteration. 
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2.6 Geometry Treatment 

Geometric modelling and grid generation is a significant portion of the effort spent 
in obtaining the flowfield about any reasonably complex geometry. A structured 
approach is utilized for this study, with the body-conformal internal grids generated 
using the elliptic techniques of Thompson, et al. [76], Thomas and Middlecoff [77], and 
Steger and Sorenson [78]. The external flow domain was discretized via algebraic and 
hyperbolic means, the topology was chosen to allow the use of these grid generation 
methods. The grids used in this investigation were generated using codes written 
by Steinbrenner, Chawner, and Fouts [79], Chan and Steger [80], and Atwood and 
Vogel [81]. A discussion of the treatment of the surface, the grid topology, and the 
grid strategy is given below. 

Surface Modelling 

The geometry used for the SOFIA configurations utilized clipped wings to emulate 
the geometry used in a specially designed experiment [16]. The use of clipped wings 
in the wind tunnel test allowed a cavity of more realistic size to be studied. The 
fuselage, wing, fairing, nacelles, and telescope geometry were obtained from CAD 
databases. Positioning errors in the database were corrected using blueprints. 

The process of generating the more complicated grids, e.g., the quiet SOFIA 
configuration with telescope (configuration 100), warrants additional comment. Ex- 
tensive wind tunnel testing [16] resulted in a hand-formed ramp and aperture, shown 
in Fig. 5a. This geometry was subsequently laser digitized [82], resulting in a data set 
of the form shown in Fig. 5b. These data, accurate to approximately 0.2 mm (0.6% 
of cavity length), were then converted into a form suitable for the surface grid using 
a standard CAD package [83] . 

Surface definition via bicubic surfaces in regions of high curvature can cause local 
oscillatory behavior [79, 81]. Along overlapping surfaces this property is manifested as 
C°, or jump discontinuities at zone boundaries. The problem is ameliorated through 
bilinear projection from one zone boundary to another [84], The distance of projection 
is typically five orders of magnitude less than a characteristic geometry length. 
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Figure 5: Geometry acquisition: (a) model, (b) laser digitized configuration, and (c) grids 
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Topology 

The overset grid topology scheme was chosen for its geometric flexibility and its ability 
to allow refinement of individual zones. The use of a different grid for each component 
of the geometry simplifies changes that will occur as the design matures. In fact, the 
geometries with cavities were built upon the clean configuration grids, providing a 
savings of many man-hours. The topology was chosen to allow rapid evaluation of 

new configurations and permit simple specification of turbulent wall and shear layer 
regions. 

Grids 

The SOFIA configuration y + values for the first grid point away from the wall were 
generally about 4.0, the farfield boundary was placed at 20 fuselage diameters, and 
the outflow 10 diameters downstream. Damping of acoustic waves at the farfield 
boundary was achieved by the use of large cells which were unable to support the 
high frequency waves. 

The clean SOFIA 747 configuration, without a cavity, was modelled using four 
grids for the half-body: one each for the fuselage, wing, wing tip, and nacelle. The grid 
point count was approximately 4 x 10 5 . The fuselage grid was refined in anticipation 
of the cavity to provide similarly sized cells in interpolation regions. 

The untreated aperture geometry, configuration 25 of the wind-tunnel test, was 
gridded by reflecting the four grid zones described above and adding two for the 
cavity. The term untreated refers to the lack of geometry modifications which can 
eliminate cavity resonance. The fuselage zonal boundaries were shifted meridionally 
to move interpolation away from the cavity region. The two additional grid zones 
consisted of an outer cavity grid surrounding the cavity region and an inner cavity 
grid which included the cavity walls and the shear layer region. The outer zone was 
utilized to isolate the cavity unsteadiness from the global solution. The total grid 
point count for this case was about 1.2 x 10 6 distributed in 10 zones. 

The treated aperture geometry, configuration 100 of the wind-tunnel test, was 
modelled by the addition of seven grid zones to the clean case: one each for isolation, 



CHAPTER 2. NUMERICAL METHOD 


38 


aperture wall, shear layer, cavity wall, telescope tub, secondary mirror, and an inner 
ramp grid. The total grid point count for this case was about 1.8 x 10 6 in 15 zones. 


Chapter 3 

Boundary Conditions 


The flow solver block implicit boundary conditions are implemented in a manner con- 
sistent with the flux split linearization described earlier. The inviscid and viscous im- 
permeable wall conditions are prescribed similarly to those given by MacCormack [85]. 
Although the following procedures are presented for a cell face which lies along a con- 
stant 77 plane, the procedure may be generalized for application to any cell boundary. 
Finally, the characteristic inflow and outflow boundaries are discussed. 

The inviscid, impermeable wall boundary condition is described for a pair of cells 
between which the surface lies, depicted in Fig. 6. In the following discussion, the cell 
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above the wall will be denoted by subscript 2, the cell below the wall by subscript 1. 
At the centroid of cell 2 the velocity is expressed as 

V 2 = u 2 \ + U2j + w 2 k 

and at the cell wall the surface normal is 

s wall = + SyJ + S z k)| wa // 

4" SyJ 4- S;k)|tva// 

where ||s|| is the vector magnitude. Hence, the velocity component normal to the wall 
is 


K»2 = || K 2 1| £*,„// 

= ( us x 4- vs y + ws z )(s x i + sj 4- s z k) wa ii 

^ 

Since, V — V t + V n} then the tangential velocity component is 


Vt 


1 2 


= V 2 -V n2 

= [tt - s x (us x 4- vs y + u;s z )]i 

+ [v - Sy(US x + VSy + tCSjJj 

4-[u> — s z (us x + vs y 4 - ws z )]k| 2 


The flow tangency condition is satisfied by V n = V t2 and V nl = -V n2 . Total energy 
and density are found from reflection as even functions through the relation RSQi = 
ER6Q 2y where E = diag[\, 1, — 1, 1, 1] and the rotation matrix is found from an 
eigenvalue problem, the eigenvectors being the rows of 


R = 


1 0 0 

0 “-(Sy 4“ Sj j s x 

0 S X Sy 

0 s z s z 

0 0 0 


0 0 

s x 0 

s z 0 

4" Sy) 0 
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A Riemann problem can then be solved at the wall to determine the flux from F j=i = 
B Q2 + B + Qi. This amounts to the wall being represented as a contact discontinuity 
by constraining the contravariant velocity to vanish. Implicitly, this results in 


6Q 1 = R~ l ER\ wall 6Q 2 


R~ l ER = 


where 

' 1 0 

0 1 — 2 ?^ 

0 —2s y s x 
0 —2 s 2 s x 

0 0 0 

The block tridiagonal system may be written as 


0 

2 S%Sy 

1 - 2 s 2 

— 2S Z Sy 


0 

-2s x s z 

~~2SyS Z 
1 - 2 5 * 
0 


0 

0 

0 

0 

1 




( 6 ) 


’ B[ C[ 


6Q1 ' 


' aq, ■ 

A' 2 B’ 2 c 2 


6Q 2 

— 

aq 2 


Now the change in flux across an arbitrary cell wall boundary is given by AE wa u = 

A+6 Q 1 + A ~ d Q^ or the sum of the changes in the flux contribution from the positive 
and negative moving waves. Substitution of Eq. ( 6 ) yields 


AE wa ii = (A + i? 1 ER + A )6Q 2 

and it can be seen that dependence upon 6 Q x has been eliminated. Hence, the block 
tridiagonal system may be represented with embedded boundary conditions as 


' B" C' 2 


6Q 2 


A Q 2 

A' B' 3 C’ 3 


6Q 3 

— 

aq 3 


where B" is the appropriately modified Jacobian. 

The viscous impermeable wall imposes additional constraints on the specification 
of the wall flux. Again utilizing a primitive variable vector V = [p, u, v, w, e] T , then 

at r 

6V X = diag[l,t u ,t vt t w ,t]6V2 = ^ 77-^2 
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for a wall face at §. In this form the toggles t u ,t v ,t w , and t are set at 1 or 1 for a slip 
or a no-slip condition, or adiabatic or isothermal wall, respectively. This may be seen 
by simply rearranging expressions of the form u wa u = 5(^1 + «2) or u waii = u i — u 2- 
Having already specified the impermeable wall conditions earlier, only the viscous 
terms at the wall are of present concern. Looking at the terms of the form 

AQ 2 = ^{-MNiSQi + MN2&Q2 H ) 


then substitution of the wall relations above leaves a term 


~ At I ,, 
AQ 2 = ^ -M 


§«!i < g 1+ Ag ! + 

dVi dQi dQi 


At 

V 2 


W2 6Q2 + ' 


which is subsequently embedded into the block structure. The dependent variables 
within cell 1 are specified according to boundary-layer theory, holding the pressure 
gradient zero normal to the wall. The remaining variables follow from fluid and 
thermodynamic relations. 

The inflow and outflow boundaries are specified according to characteristic theory 
for generality. Linearization of T^Qt+AT^ = 0 for the forward differenced inflow 
condition yields 

[/ - yA] T-HQt = -yA f-\Q^ - Q") 

where the modified eigensystem matrices are computed as 


A = diag [0, 0, 0, 0, (1 - t in )(U - c )] ||s|| 



dp t/9Q 

dTrjdQ 

dv/dQ 

dw/dQ 

(1 — t, n )^5 +ti n du/dQ 


Here t in is zero or unity for subsonic or supersonic inflow. The specified variables are 
chosen such that a unique set of flow quantities are given at the entrance. 
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The outflow condition is specified in a like manner; however, there are now at least 
four characteristics linking the domain with the boundary. Backward differencing 

about i + 1 and specification of static pressure at the exit results in the modified 
matrices 

A = diag [ 0 , U,U,U 4- c, t^U - c)] ||s|| 

h 

h 

f- 1 = h 

u 

. (1 - tout)dp/dQ + t out l b 

In the above development, the eigensystem is evaluated at the boundary face in 
question, maintaining consistency with the interior treatment. 

The strongly unsteady blast-wave problems investigated here revealed that the 
use of block implicit boundary conditions resulted in significantly enhanced conver- 
gence. This beneficial effect is caused by the faster signal propagation arising from 
the incorporation of the boundary conditions within the linear system. However, for 
the cavity flows explicit boundary condition implementations were used for coding 
simplicity. Comparison of the computed results with experiment show satisfactory 
resolution of the moderate unsteadiness present in transonic cavity flows. 

Several of the cavity cases used characteristic boundary conditions holding mass 
flow, total enthalpy, and flow angle constant. Subsonic inflow conditions, for example, 
are related to the interior of the flow by the u — c characteristic: 



which may be rewritten as 

p t — pcu t -f c(pu t 4- up t ) = — (u — c)\p x — pcu x ] 4- c(pu) t 

where the mass flow is fixed for these computations, giving (pu) t = 0. Discretization 
and rearrangement yields 

j — (u — c)^- 

dp = ~~ 3e T [P 2 ~ Pi ~ pc{u 2 - «i)] 

dp ■+■ cu 
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Implementation of the boundary conditions, unless otherwise noted, are as follows: 
viscous impermeable wall conditions are no-slip, zero normal pressure gradient, and 
adiabatic; information transfer across overset mesh boundaries is implemented using 
non- conservative trilinear interpolation of Q. TVeatment of the farfield boundaries is 
case dependent and is noted in the results section. 



Chapter 4 

Geometrical Aero- Optics 


The objective of the numerical simulation of the flow about the SOFIA airborne ob- 
servatory is to design a safe configuration which will have the least detrimental effects 
upon the optics. Towards this goal, the following transonic cavity flow problems were 
divided into three sections. First, the unsteady interaction of the external flow with 
the cavity requires time-dependent solutions to the Reynolds-averaged Navier-Stokes 
equations. Second, the shear layer growth rate is strongly dependent on turbulence 
effects which must be modelled due to the grid coarseness. The final portion of the 
problem is the application of the optical model to the unsteady density field in the 
shear layer to determine seeing quality. 

The variation of the speed of light through gases is primarily a function of the 
density field. This fact has been extensively used to benefit the study of fluid physics, 
as exemplified by use of schlieren, shadowgraph, and interferometry techniques. How- 
ever, the objective of the present effort is to quantify the wavefront distortion of a 
beam of light propagating through the shear layer. This distortion is computed using 
the history of the density variations within the shear layer to predict fluctuations in 
the optical path length via geometric optics. This will in turn allow prediction of the 
telescope resolution limits due to seeing and thus contribute to the telescope design 
specifications. 

The geometric optics model developed here assumes that the impact of the fluid 
density on the optical field may be computed by casting light rays through a field 
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Figure 7: Geometric optics: (a) partitioning of hexahedrons into tetrahedrons, and (b) 
intersection and refraction procedure 


discretized into tetrahedrons. Diffraction effects, which become important when the 
wavelength approaches the scale size, are neglected. The simplifications afforded by 
the use of planar facets and piecewise continuous media are utilized by tesselating, 
or partitioning, each hexahedron of the flowfield into five tetrahedra as shown in 
Fig. 7. It should be noted that other tesselations were found to be more robust for 
thin warped cells [86], however for the shear layer grids used here the five tetrahedra 
decomposition is well-behaved. 

Application of the geometric optics code to two preliminary test cases was under- 
taken to determine sensitivity of the optics code to the above non-unique tesselation. 
The parallel emergence of the rays after propagation through a plate and a prism of 
index of refraction n = 2.4 suggests that the results are relatively insensitive to the 
method of tesselation. 

The problem can now be divided into three steps: 1) propagation of the ray, 2) 
intersection of the ray with a facet, and 3) refraction. Solution for the point contained 
in both the facet, p, and the ray, q;, expressed parametrically as 

q,(t) = d + et ; p(u, w) = a + bu + cw 

results in the intersection in parametric coordinates, which is found from the dot 
product of the normal with the ray and the surface: 

(b x c) • p = (b x c) • q, 

_ (b x c) • a — (b x c) • d 
(b x c) ■ e 
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(c x e) • d — (c x e) • a 
(c x e) • b 

(b x e) • d - (b x e) • a 
(b x e) • c 

where the vector coefficients are found from boundary conditions: 

a = p(0,0) 

b = p(l,0) — p(0, 0) 

c = p(0, 1) — p(0,0) 

d = q,(0) 

e = q.(i)-q,(o) 

Specification of the light ray origin and a direction initializes the problem. Following 
the search for the initial hexahedral cell in which the ray originates, the tetrahedron 
within this hexahedron must be computed. First, the shortest intersection distance 
of the ray with the 16 planes which compose the hexahedron is computed. Then 
the dot product of the ray with the fourth vertex of the closest plane determines the 
origin tetrahedron. Subsequent intersection and refraction processes are a marching 
procedure. The optical path length {OPE) is found from 

OPL = J n(s)ds « ^ n, As, 

j 

where n(s ) is the index of refraction as a function of position along the ray, s. The 

variation of the OP L over the aperture gives a measure of the wavefront error caused 
by the shear layer. 

The refraction process is determined according to Snell’s law as shown in Fig. 7, 
where the planar interface, p(u, w), separating the media and the incident light ray, 
q,(t) [87, 81] are depicted. Generalization to three-dimensions is accomplished by 
rotation to the osculating plane, which includes the surface normal and both the 
incident and refracted rays. In this osculating plane, a rotated local coordinate system 
is defined: 


q. = |q<Jfi + |q.,|t 
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where 

|q.„l = |q.|costf, and |q, t | = |q,| sin 6>, 

- _ q, — n cos 0, 
sin $i 

Application of Snell’s law n\ sin 0, = «2 sin 0 r where n = 1 4- results in an 

expression for the refracted ray: 

q r = n cos 9 r + 1 sin 6 r 

The local index of refraction, nj, is found by arithmetically averaging the densities 
at the four vertices of the tetrahedron, where only one vertex changes as the ray 
propagates to a neighbor tetrahedron. The Gladstone-Dale constant, P, is a function 
of the media and of the wavelength. Using air as the media and a wavelength of 
A = \ D = 5893A, then P = 2.92 x 10 -4 . The Cauchy formula can also be used: 

(3 = [2875.66+ 13.412/(A 2 x 10~ 8 ) + 0.3777/(A 4 x 10" 16 )] x 10 -7 

The values of P used in the present computations were chosen to match those used in 
the reduction of the experimental data. The wavelengths of interest for SOFIA range 
from the near infrared, l/^m, to the microwave, 1 mm, where optical distortion can 
be seen to be more severe for shorter wavelengths. 

Finally, to obtain a measure of the loss in irradiance due to the fluctuating density 
field, the OPL for vacuum conditions is subtracted from the OPL through the gas 
to yield the optical path difference ( OPD ). The value of < OPD' > is computed 
using a sequence of OPD ' s at a fixed station. Using the root-mean-square wavefront 
distortion < OPD ' >, the phase distortion ($) is found from $ = ^22. The Strehl 
ratio, given by 



is a measure of the peak intensity to which a beam can be focused. The computational 
expense of the procedure outlined above is currently 250/xs/hexahedron/ray on a 
single processor of the Numerical Aerodynamic Simulator Cray 2. 

In these studies, the effect of fluids upon the optical field is determined through 
prismatic modelling of the density field. Integration of the equations of motion 
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through a trilinearly varying density field could also be implemented. The use of a 
six-tetrahedron decomposition would eliminate the present requirement of a checker- 
board cell arrangement to prevent gaps between hexahedrons. The modelling of the 
wave-like nature of light could also be implemented via the parabolized Helmholtz 
equation if diffraction effects were deemed significant [88]. 



Chapter 5 

Results and Discussion 


5.1 Blast- Wave Results 

The methods introduced in the previous sections are applied to test cases which 
demonstrate the capabilities of the algorithm. The viscous term treatment in a low 
Mach number regime is shown by the Couette flow problems, which are compared 
to Similarity solutions and previously obtained numerical results. Demonstration of 
the inviscid term treatment is shown by capturing of transient discontinuities for a 
shock tunnel start-up problem. The three-dimensional results are compared with an 
experimental study of a hemicylinder mounted in a shock tube. 

5.1.1 Couette Flow 

The Couette flow problem is used to compare the present methods against the method 
of Beam and Warming [7] and the similarity solution as given by Schlichting [89]. The 
results for the two upwind methods fall virtually on top of each other, and n indicates 
the time step. The solutions shown in Fig. 8 were obtained using quiescent initial 
conditions and viscous boundary conditions with no-slip adiabatic walls. Both of 
these cases were implemented in the thin-layer form at a Reynolds number of 6.4, 
based on the distance between the plates, equal to 10 -5 feet. During the course of 
these solutions, slightly more than an order of magnitude drop in \Sp\ max per two 
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subiterations was observed in the (3 x 10) cell domain. The Courant number used 



Figure 8: Couette flow case: (a) impulsively started and (b) oscillating plate 

for the oscillating plate calculation was approximately 10, indicating the viability of 
these types of unsteady computations at Courant numbers greater than unity. The 
Courant number was computed using CFL = fma X [\(U y V ,W)\ + c]||s|| over each 
cell in the domain. Identical results were obtained using both the two- and three- 
dimensional implementations in all directional permutations. Results reveal slightly 
steeper gradients than that of the conventionally differenced scheme or the analytic 
solution, a possible consequence of the handling of the boundary conditions or the 
viscous term treatment. In addition, this case was found to be insensitive to the 

choice of the higher-order flux correction terms, possibly because of the dominance 
of diffusive effects. 

5.1.2 Shock Tunnel Start-up Problem 

The third test case evaluated the inviscid term treatment through the simulation of 
the transient starting process of a planar shock tunnel. The (300 x 60) cell domain 
is shown in Fig. 9. The solution of the Euler equations is presented in Fig. 10 as 
a comparison of experimental and numerical shadowgraph images, the former due 
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to Amann [90]. The computed shadowgraph function is proportional to V 2 p and 
thus acts as an amplifier of the density gradient across contacts and shocks, for in- 
stance. This solution was obtained using Roe flux difference splitting with <j> = 1/3, 
the upwind biased flux evaluation. The Steger- Warming flux evaluation with Roe 
averaging was found to be moderately less stable, but no significant differences in the 
results were found for this case. The maximum compression parameter was used, and 
the entropy-fix parameter used in Harten’s formulation [9] was set to 0.15. Discon- 
certingly nonphysical solutions were produced for smaller entropy-fix levels, possibly 
associated with an entropy- violating condition. The problem was initialized with a 
moving shock propagating to the right at a Mach number of 2.97, while the boundary 
conditions were specified as impermeable inviscid along the walls and fixed for the 
inlet and exit. For the maximum Courant number of four used here, four subiter- 
ations were chosen per time step based on a subjective judgment of discontinuity 
sharpness. The \6p\ max was observed to drop approximately an order of magnitude 
over the course of these four subiterations. 

Physically, this nozzle starting process generates a high enthalpy reservoir of more 
than 50 times the initial pressure, while the density increases by 11 times the initial 



CHAPTER 5. RESULTS AND DISCUSSION 


53 



Experiment, Amann (Ref. 90) Euler 


(a) Primary shock reflected into reservoir 

(b) Swallowed primary shock 

(c) Rearward facing shock being swept downstream 

(d) Reflected shock system 

(e) Mach line generated from C 2 discontinuity 


Figure 10: Shock tunnel case: shadowgraph comparison for a Mach 3 planar shock 
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state. This reservoir provides the energy necessary to generate high Mach flows 
downstream of the diverging nozzle region for short durations. The ensuing reflections 
of the shock with the nozzle wall reveals the complexities of the shock-shock and 
shock-contact interaction. In particular, it can be seen that the development of the 
rearward facing shock, which is directed upstream while being swept downstream, 
is resolved. At later times, the finer scale fluid motion between the primary and 
rearward facing shocks is, for the most part, lost because of grid coarseness and 
attendant numerical dissipation. However, increasingly fine structures are captured 
as the grid is refined. 

5.1.3 Shock Tube Blockage Study 

The viability of the technique in three-dimensions is shown by the final test case. 
These results are intended to replicate the conditions in an experimental study of 
a blast- wave encounter with a hemicylinder target in a shock tube by Kingery and 
Bulmash [91]. The experimental test configuration and pressure transducer locations 
are shown in Fig. 11. In order to estimate the costs and benefits of inviscid versus 
viscous simulations, the flow about this geometry was computed using both the Euler 
and Navier-Stokes equations. However, the expense of these three-dimensional simu- 
lations permitted the use of only one of the inviscid flux evaluation methods; the Roe 
flux difference splitting was chosen. 

The simulation was initialized as a translating planar shock before diffraction over 
the cylinder began. Initial conditions were specified as: 

M, = 1.518, Re/m = = 23.3 x 10 6 /m 

^ = 288.17 K , poo = 101.3 x 10 3 N/m 2 

Boundary conditions are specified for the viscous, single zone computation as shown in 
Fig. 12. In the shock-tube direction, the f -direction, extrapolation is used. This non- 
physical extrapolation is adequate for the duration of the early interaction. However, 
solutions at larger times are suspect, where times after the shocks have propagated 
through the boundaries are defined as large. Additionally, the use of an advancing 
front boundary is enabled because of a priori knowledge of the grid structure and the 
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Figure 11: Hemicylinder case: experimental configuration and pressure transducer loca- 
tions 
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Figure 12: Hemicylinder case: central region of the 78 x 50 x 25 cell (a) inviscid and (b) 
viscous grids 


primary shock speed. This simple time-dependent boundary reduces the computation 
time by using the fact that nothing occurs ahead of the blast- wave. In the 77-direction, 
the lower boundary defines the surface geometry of the hemicylinder, and hence is 
specified as a no-slip isothermal wall. The top of the domain in the 77-direction, corre- 
sponding to the inner radius of the shock tube, is specified as an inviscid wall, based 
on the assumption that the viscous effects on this surface have negligible influence on 
the results. Finally, the (-direction boundaries are treated using the viscous condition 
along the floor of the tube. Symmetry conditions are used along the plane running 
along the longitudinal axis of the cylinder and normal to the floor. To simulate ex- 
perimental conditions, the wall temperature was set equal to the temperature of the 
quiescent flow prior to primary shock arrival. The viscous grid has normal spacing 
of approximately 10 — ^ meters at the viscous walls. The Euler computation used the 
inviscid boundary conditions previously discussed where appropriate. For this Euler 
grid, since the areas of the faces corresponding to the geometric axis singularity are 
zero then F ■ s is also zero. Compensation for the round-off error inherent in the grid 
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was implemented by eliminating those face areas which fell below a specified toler- 
ance. The grids and boundary conditions for these cases are partly shown in Figs. 12 
and 13. 


iHiiiliniSM 



Figure 13: Hemicylinder case: symmetry plane of the central part of the viscous grid 


The inviscid computation used <t> = 1/3, /? = 4, Harten’s entropy fix parameter of 
10 4 , Roe flux difference splitting, one subiteration per second-order accurate time 
step, and a Courant number of 15. The solution was obtained in 1500 time steps 
without any change of parameters. 

The viscous computation used the same flux evaluation as above with the addition 
of the second-order accurate full viscous terms. Because of the viscous spacing, the 
Courant number utilized was 10 4 , allowing the solution to be obtained in 6800 time 
steps with no subiterations. In contrast to the inviscid simulation, the advancing 
front boundary condition was utilized in this case. 

Results, given in Figs. 14 through 17, show that the primary shock is captured 
over two to three cells, the large physical thickness obtained is an artifact of the 
coarse grid used. Adaptive gridding methods would help maintain sharp shocks, but 
the anticipated expense of these methods precluded their use here. Figures 14 and 15 
show comparisons of the numerical and experimental pressure histories. It is seen that 
the peak overpressure is underpredicted by 10%, possibly owing to the coarseness of 
the grid which in turn thickens the shock. These computed surface pressures were 
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extracted from the domain through the use of a Newton search in three-space for 
the cell in which the given (x,y, z) probe coordinate fell [92], Subsequent trilinear 
interpolation over the cell, where the uniform parametric coordinates ( u , v, w ) are 
determined from the positions of the vertices of the hexahedral cell, allows the pressure 
to be computed. Inherent in this first-order approximation lies the assumption that 
over a discrete cell the variation of pressure is linear in space. 

Figures 16 and 17 show portions of the viscous simulation at selected times. Phys- 
ically, the interaction process begins with the normal impact of the incident shock 
with the front face of the hemicylinder. At this time, peak overpressures of six times, 
and densities of four times that of the quiescent state are generated along this for- 
ward face. As the shock diffracts over the sharp corner of the target, a separation 
bubble forms, which eventually envelops a large portion of the circumferential face of 
the body. This vortical motion is depicted in Fig. 17 by instantaneous streamlines. 
A supersonic pocket is generated as the air negotiates the sharp corner as it rushes 
from the stagnation region left in the wake of the upstream propagating reflected 
shock. The next significant event occurs as the shock diffracts over the rearward face, 
shedding a strong vortex sheet while an expansion wave propagates away in a pattern 
which grows with time. The diffracted shock then impacts the floor of the shock tube, 
reflecting it upwards, while the shock which diffracted over the circumferential face 
reflects inwards from the outer walls of the tube. A simplified sketch of the interac- 
tion process is shown in Fig. 18. The subsequent diffractions and reflections result in 
the interaction of shocks, expansion fans, vortices, and developing boundary layers. 
From experimental evidence, this gross unsteadiness does not dissipate for more than 
15 milliseconds after the interaction event begins. However, the primary shock passes 
from the test section 5 milliseconds after the initial target interaction; therefore, the 
computation is stopped at that time. 

The effects of the viscous terms are seen by comparing the pressure histories in 
Figs. 14 and 15. While the pressures along the upstream face are largely unchanged, 
the circumferential and downstream faces are significantly affected by viscosity. The 
large separation along these faces causes low pressure regions due to this vortical 
motion. This phenomenon is more accurately captured in the viscous simulation, as 
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Figure 14: Hemicylinder case: experimental [91] and inviscid computation pressure histo- 
ries 
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Figure 15: Hemicylinder case: experimental [91] and viscous computation pressure histo- 


ries 






Figure 16: Viscous hemicylindercase: symmetry plane pressure [atm.] and Mach contours 
at (a) t=l.lms, (b) t=1.7ms, and (c) t=6.2ms 
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Figure 17: Viscous hemicylinder case: instantaneous streamlines at (a) t=l.lms, (b) 
t=1.7ms, and (c) t=6.2ms 
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Figure 18: Hemicylinder case: schematic of shock interaction 


may be seen by inspection of the pressure histories at probe 11. Differences between 
the experimental and the present results may be due to poor capturing of the vortex 
strength owing to grid coarseness. However, the higher-order behavior of the method 
used here attempts to reduce the need for finely spaced meshes. In addition, the 
occurrence of deformation of the shroud wall is thought to be a possible event during 
the experiment, and could adversely affect the comparison between experiment and 
computation [93]. 

A limited cost /benefit study of the Euler versus Navier-Stokes equations was also 
performed for the hemicylinder case. For approximately 5.5 ms of flow history on a 
(78 x 50 x 25) cell grid, the Euler computation consumed 7.6 processor hours, while 
the viscous simulation required 18.2 hours. From these results, the somewhat more 
accurate solution given by the Navier-Stokes simulation may be worthwhile. This 
is particularly true if the flowfield behavior after the direct interaction the primary 
shock with the target is important. 
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5.2 Cavity Flow Results 

Validation of the diagonalized code was accomplished by evaluation of two- and three- 
dimensional cases related to the transonic aero- window problem. Numerical results 
for free shear layer and rectangular two-dimensional cavity flows were compared with 
analytic and experimental data to evaluate the capability of capturing the fundamen- 
tal physics. Analysis of the SOFIA configuration simulations, including evaluation of 
optical distortion is also presented. 

5.2.1 Free Shear Layer 

A series of numerical experiments was performed using a two-dimensional shear layer 
as the test case. Sensitivities of mean and time-varying quantities to changes in time 
step size, fourth-order dissipation levels, and grid refinement were determined. Addi- 
tionally, partial validation of the algebraic turbulent shear layer model was determined 
through comparison with similarity solutions and experimental data. 

The computational domain for this case includes a two inch long splitter plate 
embedded in a channel, with initial conditions specified as a discontinuous step at 
the channel centerline. Shown to scale in Fig. 19, the channel extends 30 inches 
downstream of the splitter plate trailing edge, and five inches above and below the 
plate. Inviscid walls were specified for three inches upstream of the viscous splitter 
plate and for the channel walls. The inflow and outflow conditions were implemented 
using one-dimensional characteristic relations holding mass flow, total enthalpy, and 
flow angle fixed at the inlets, and fixing pressure at the exit plane. The boundary 
layers on the splitter plate and the shear layer were turbulent. Reynolds number 
based on the mean velocity of the streams and the length of the splitter plate was 
6.7 x 10 5 . 

The results for three grid refinement levels are shown in Fig. 19 along with 
Gortler’s similarity solution. The velocity profiles are taken 10 inches downstream of 
the trailing edge of the plate. The solution can be seen to become grid independent 
when the grid becomes finer than approximately 20 points across the layer. The Mach 
number ratio for this case was 0.2/0. 8 and a = 20.7. Eddy viscosity was observed to 
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Figure 19: Shear layer case: (a) velocity profiles of differing grid resolution compared to 
similarity and (b) variation of spread rate with velocity parameter 

grow linearly in accordance with the Clauser formulation. 

Numerical experiments to determine the dependence of < p' > on the level of 
fourth-order dissipation showed that a change in fourth-order smoothing from 0.01 
to 0.05 caused a change of less than 1% in sound pressure level. 

The velocity ratio across the SOFIA cavity shear layer will vary with streamwise 
location. Hence, comparison of the variation of spread rate with velocity ratio is 
shown in Fig. 19b. Three velocity ratios are shown for low Mach numbers with about 
20 points maintained across the layer for all cases. Although the computed spreading 
rates are within the bounds of the experimental data [59], the data point at r = A 
falls below the trend because of the limited entrainment afforded by the inviscid side 
wall treatment. 

5.2.2 Two-Dimensional Rectangular Cavity 

The objective of this two-dimensional cavity case was to demonstrate the prediction 
of self-induced cavity resonance. Validation data are provided by comparison against 
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Rossiter’s experiment [15]. Sensitivity of the solution to topology, second-order dis- 
sipation, and turbulence model effects were determined. 

The cavity geometry and grid topology are shown in Fig. 20, where the grid has 
been coarsened for clarity. The test conditions were set as: 

Moo = 0.9, Re i = 1.47 x 10 6 , 1 = 8 in. 

Poo = 0.40 kg/m 3 , = 2.9 x 10 4 N/m 2 

The ratio of cavity length by depth ( L/D ) was 2 for this model. The inflow boundary 
was placed 7.5 L upstream of the cavity leading edge, the outflow boundary 4.5 
L downstream of the cavity trailing edge. The inflow and outflow conditions were 
specified as for the free shear layer cases, and an inviscid wall was placed 5 L above 
the cavity. 


Global Zono 




Figure 20: 2-d cavity: topology and grids 


Figure 21 depicts instantaneous Mach number and pressure contours obtained 
during the computation in which the time step size was At = 1.97 ps. Inspection 
of the contours across zonal boundaries indicates that the interpolation process is 
well-behaved for this unsteady flow. The Mach number contours show an instant 
of the time-oscillatory shear layer behavior which is prevalent for these rectangu- 
lar geometries. The pressure contours verify the feedback mechanism postulated by 
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Figure 21: 2-d cavity: instantaneous Mach number and pressure [kPa] contours 
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Rossiter [15], and shown graphically in Fig. 22. Briefly, the cycle begins with the 
propagation of a wave from the aft wall of the cavity to the foiward face. Wave 
reflection from the forward wall causes the shear layer to bow outwards, shedding 
vorticity. The deflected shear layer convects downstream and induces another cycle. 

The origin of this physical model can be linked to the edge-tone phenomena, where 
a thin planar jet interacts with a wedge. The frequencies at which this feedback is 
reinforced is determined by ambient temperature and Mach number and was first 
quantified by Rossiter. Derivation of this model begins with the assumption that the 
frequency of the vortex shedding is equal to the cavity acoustic field, / = -£*• = yf • 
The vortical and acoustic field relationships are linked by 

m„A„ — L + 7 „A v + Kuoot' 

L = m a A„ + c T t' 

from which 

, _ Uqq (m-7) m 

L (£ + 7 ?) 

where the phase lag factor, 7 = 7. = 0.25, and the normalized convection velocity 
of the perturbations, K = 0.66, are empirically determined constants dependent on 
the geometry and ambient conditions. The integer stage number is given by m = 
m a + m v . Use of the a Mach number scaled by the cavity speed of sound, determined 
by the recovery temperature, offers improved correlation with experiment. The model 
described here is idealized: the shear layer perturbations may be manifested as a 
sinuous motion as opposed to discrete roller vortices depicted in Fig. 22. 

The pressure histories along the cavity walls are depicted in Fig. 23 along with 
the comparison of Rossiter’s data to present results in power spectral density ( PSD ) 
form. The comparison against experiment shows agreement in frequency at the peak 
magnitudes. Magnitudes are higher for the present case by about 2 dB, however 
for this experiment and other numerical work this has been observed as an effect of 
cavity width. In these studies [15, 36], the sound level was found to be inversely 
related to L/W\ the cavity for Rossiter’s experimental work was of L/W = 2. Also 
shown in Fig. 23 are the first four stages, rrij, predicted by Eq. 7 using both his 
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Figure 22: Rossiter’s feedback model 


vortex convection ratio of A — 0.66 and a K — 0.56 inferred from the treated cavity 
computation discussed below. The PSD for this case was computed using 8192 time 
samples, no zero-padding, and a square window. 

Variation of the second-order dissipation, e 2 , from 0.5 to 0.3 caused no discern- 
able change in the pressure histories. Additionally, in order to test a hypothesis of 
a limited domain of unsteadiness, an isolation zone was implemented as shown in 
Fifi- 20. The flow outside the zones of interest was frozen, resulting in a decrease 
of the cavity sound pressure levels by 2%. A final comparison between experimental 
data and numerical results is provided in Fig. 24, where the variation of the mean 
and oscillatory pressures along the cavity walls is shown. The < C > is computed 
following Rossiter’s reduction of experimental pressure histories 


[A 



e = 0.155 


where P and R indicate peak and background pressure levels. The computational 
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Figure 23: 2-d cavity: pressure history and power spectra comparison 
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Figure 24: 2-d cavity: variation of mean and oscillatory pressures 


results were ^educed from assuming p e |* = p and p ( \ P = p. Use of Reynolds averaging, 
specifically fg — / g, gives a result consistent with the experimental reduction 



Given the difference in spatial dimensions and turbulence modelling uncertainties, 
the trends shown in Fig. 24 appear reasonable for a flow of this complexity. 

Finally, to give a qualitative comparison of numerical and experimental [17] re- 
sults, schlieren images are shown in Fig. 25. Despite Reynolds and Mach number 
mismatches and grid coarseness away from the cavity, the observed and computed 
acoustic radiation patterns are similar. 
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Figure 25: Comparison of experimental and numerical schlieren images with knife edge 
vertical 

5.2.3 Two-Dimensional Treated Cavity 

The effect of cavity geometry, particularly modification of the shear layer attachment 
region, is known to possess potential quieting capabilities [42]. The Army Airborne 
Optical Adjunct (AOA), shown in Fig. 26, flight tested several passive and active qui- 
eting methods [94]. The purpose of the present numerical simulations is to determine 
if optical quieting methods, particularly aft ramp treatment and lip-blowing, could 
be accurately simulated. Tangential lip-blowing at the upstream edge of the aperture 
may provide quieting by replenishing the mass entrained by the shear layer from the 
cavity. The quieting provided by aft ramp treatment at the shear layer impingement 
region is discussed below. 

The grid cell size was specified at 0.83" in the streamwise direction, chosen so 
that frequencies up to approximately 400 Hz would be resolved without significant 
numerical dissipation effects. A time step of 44/xs was fixed so that CFL ~ 1 in 
the streamwise direction within the shear layer. The numerical test conditions were 
matched to flight data: 

Moo = 0.77, Re L = 5.00 x 10 6 , 

Poo = 0.262 kg/m 3 , Poo = 1-63 x 10 4 N/m 2 


L — 47 in. 
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Figure 26: U.S. Army Airborne Optical Adjunct [94] 


The initial conditions used the assumption of isentropic recovery to obtain the cav- 
ity temperature while maintaining constant pressure across the aperture. Boundary 
conditions were of the one-dimensional characteristic form: constant mass flow, total 
temperature, and flow angle inflow; constant pressure outflow; and an inviscid wall 
6.4L from the cavity. A characteristic inflow condition was also used for the lip- 
blowing boundary, the flow rate computed using flight data and assuming isentropic 
compression of the ram air utilized in the aircraft. The 100% lip-blowing rate case 
corresponded to a m = 0.42(pu) oo . For the discussion below, computed high and 
low lip-blowing rates refer to 100% and 1% of this mass flow rate. The coarsened 
near-field grids are shown in Fig. 27. 

Comparison of flight data with this planar numerical study is justified by the flight 
test effort toward establishing two-dimensional flow across the apertures of the AO A. 
Use of flow cones and a shear layer rake verified, for the most part, the success of this 
effort. Although the cupola which allows for the cavities is of hemicylindrical form, the 
center of rotation is through the center of the cavity volume. Rather than simulate an 
axisymmetric cavity with an erroneous radius of curvature, the cavity was modelled 
with no spanwise variation in the flow. However, there is a dimensional effect in that 
the mass removed from the cavity by the shear layer entrainment process can only 
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Figure 27: 2-d treated cavity: near field grids 

be replenished at the impingement region. This is in contrast to the mass addition 
mechanism present in three-dimensions, which also includes mass replenishment via 
spanwise structures such as streamwise vortices. 

The mechanism by which an aft ramp reduces the cavity feedback was explained 
by Heller and Bliss [24]. Assuming two-dimensional incompressible flow, the region 
immediately surrounding the stagnation point can be treated using a streamfunction 

approach: 

= axy + ]rby 2 

M . 

ti = — = ax + by 

ay 

v = -«=- = -ay 

OX 

and the resultant field is shown in Fig. 28. 

Physically, this result can be explained from a force balance normal to a streamline 

approaching the stagnation point. About the stagnation point, the velocity gradient 
across the impinging shear layer creates a pressure gradient. However, there is a 
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Figure 28: Streamfu notion near a stagnation region 


counteracting pressure gradient, pq 2 /r, due to the differing radii of curvature above 
and below the dividing streamline. 

For a rectangular cavity, the extreme of normal impingement of the shear layer 
onto the aft bulkhead causes further deflection into the cavity. Mass ingestion into 
the cavity causes increased pressure, deflecting the shear layer outwards. With the 
shear region now outwardly deflected, mass expulsion from the cavity reduces the 
cavity pressure, inducing another cycle. Therefore, between the extremes of a normal 
or tangential impingement of the shear layer, a balance of forces may be found. Use 
of a ramp instead of a convex surface at the reattachment region prevents shear layer 
perturbations from inducing instabilities of the type seen in rectangular cutouts. The 
length of the ramp must be large enough to accommodate the magnitude of the 
transverse shear layer excursions expected during operation. 

It is hypothesized that the use of a modestly concave surface at the impingement 
point may provide increased quieting. Improved stability may be provided by the 
following reasoning. As the shear layer is perturbed downwards, the streamlines below 







CHAPTER 5. RESULTS AND DISCUSSION 


76 


the impingement point will have a smaller radius of curvature while the streamlines 
above the shear layer will be less curved. The resultant pressure gradient will drive 
the shear layer upwards. Conversely, as the shear layer deflects upwards the steeper 
tangent angle creates a smaller radius of curvature above, forcing the shear layer 
downwards towards the nominal stagnation point location. 

Computed and flight mean Mach number profiles are compared in Fig. 29 for 
two lip-blowing rates. The quantity <p indicates the angle from the cupola crest at 
which the data was measured. Figure 29 also shows Mach number contours for the 
two lip-blowing rates above each set of profiles. The Mach number contours are 
instantaneous while the profiles were averaged over 2000 time steps. The difference 
between experiment and computational results on the lower edge of the shear layer 
(r < -2") may be due to blockage in the cavity of the aircraft. This blockage was not 
computationally represented due to the lack of a complete geometry description. The 
difference at the upper aft portion of the shear layer (r > 2") appears to be due to 
blockage in the computational model. Overall, the maximum vorticity as a function 
of x-station is in agreement for both cases. 

Comparisons of power spectra at the aft ramp are shown in Figs. 30 and 31. The 
computed spectra can be seen to be quantitatively and even qualitatively different 
from flight data. The computed result lies more than 15 dB below the data, and 
a peak in the low lip-blowing rate spectra is clearly computed, but is not seen in 
the flight data. The power spectra were computed using 4096 points and a square 
window. Figure 32 shows mean and fluctuating quantities along the cavity walls, with 
the available data allowing comparison only along the aft ramp. The computed trend 
in < C{, > is in agreement with measured data, but a large discrepancy in magnitude 
is evident. The large spanwise variation in measured < C' > indicates the existence 
of three-dimensional effects or experimental errors. 

It has been noted from experimental evidence [95] that the frequency of large struc- 
tures in shear layers is independent of axial station and occurs at Strouhal number 
of st = = 0.024 ± 0.003, where 6 is the local shear layer momentum thickness 

and / denotes frequency. This phenomena is corroborated by the reduction of other 
researchers’ data [48, 49, 51, 96, 97], whose results range from about St = 0.02 to 
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Figure 29. 2-d treated cavity: instantaneous Mach number contours and mean profiles at 
(a) low and (b) high lip-blowing rate 
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0.03 for incompressible shear layers. 

Gortler’s solution can be used to obtain 6 = 0.036y^x for ao = 11.0, which 
compares favourably to the empirically determined correlation [98] of 0 = 0.034 
Using this relationship along with a compressibility correction [99], the computed 
peak in the AOA solution at 340 Hz corresponds to a Strouhal number of 0.032. For 
comparison, the peak at approximately 1800 Hz in the quieted SOFIA case corre- 
sponds to a St = 0.030, as will be shown in section 5.2.7. It is interesting to note 
that by using Rossiter’s formula, Eq. 7, the frequencies obtained for m = 4 and 5 are 
285 and 360 Hz, respectively. From Fig. 34 it can be observed that m v = 4, implying 
that m a = 1. 

Based upon these observations, it is hypothesized that large scale shear layer struc- 
tures are beginning to be resolved. However, the lack of empirical support from the 
flight data pressure power spectra is at odds with this supposition. The comparison 
is further clouded by the reasonable comparison in < p' > for the low lip-blowing rate 
shown in Fig. 33. The discrepancy may be caused by three-dimensional effects, angle 
of attack sensitivity, or geometry simplifications. 

In this author’s opinion, three-dimensional effects are the most plausible expla- 
nation for the discrepancy. Rockwell [100] noted that for sufficiently large Reynolds 
numbers three-dimensionality reduces coherence in the shear layer. This implies that 
assumption of two-dimensionality for small flow oscillations may be suspect. The evo- 
lution of streamwise-oriented vorticity interacting with the primary vortices would act 
to spread peaks in the reattachment ramp pressure spectra. 

As a final note, Fig. 30 also depicts data, the ordinate scaled by and the 

abcissa by obtained from an AOA wind tunnel test [101], The data can 

be seen to agree more closely with the computed results than with flight data, and a 
small peak exists where expected according to the above analysis. 

5.2.4 2-D Treated Cavity: Aero-Optical Effects 

Computation of aero-optical parameters requires the use of the unsteady density field. 
Figure 33 shows the computed and experimental profiles of the root mean square of 
the density fluctuations. Levels of < p' > were computed over a time segment of 
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Figure 30: 2-d treated cavity: power spectra, low lip-blowing rate 
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about 90 ms in increments of 0.44 ms. Using the elapsed time for a particle to 
convect across the aperture at the mean shear layer speed as a characteristic time, 
T e = jfc-, then the optical computation was taken for about nine T c . In Fig. 33, 
7 is the rake angle from horizontal, with the axis of rotation offset from the cupola 
centerline. Determination of the systematic error band on the experimental result is 
discussed below. The low lip-blowing rate result underpredicts the magnitude of the 
peak in ^ however the peak location is in fair agreement. The computed results for 

* Poo 

the high lip-blowing rate compare poorly to experiment, possibly due to inadequate 
grid resolution or the increased flow complexity of the merging shear layers. This 
type of active control is presently not a design option for SOFIA, therefore further 
effort toward improvement of the high lip-blowing case was not warranted. 



Figure 33: 2-d treated cavity: profiles with (a) low and (b) high lip-blowing 


Further investigation of the low-blowing rate case revealed the presence of large 
convecting structures associated with the shear layer. Figure 34 shows a contour plot 
of jL depicting the growth and propagation of these sinuous motions in the shear 
layer. Also depicted in Fig. 34 is a schematic of the optical model, with the initial 
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and final stations of the optical path integration are given by r 0 and rj. The large 
structures, associated with a 0.03^ vertical velocity component, are the primary 
contributors of the computed density fluctuations of the shear layer. The speed of 
the waves, as determined from Fig. 35, is 0.56uoo, below the value of 0.66ttoo inferred 
by Rossiter for rectangular cutouts, yet above the 0.51«oo determined analytically by 
Roscoe and Hankey [102]. 

Chew and Christiansen [50] and Tsai and Christiansen [51], utilizing results from 
computation and experiment, deduced that a free shear layer model of a sinusoidal 
phase delay growing in x would produce results similar to those observed. Figure 35 
displays behavior of a similar nature for the aero-window problem modelled here. 

Comparisons of integrated aero-optical quantities, shown in Fig. 36, reveal slight 
overprediction for the low lip-blowing case and, given the < pi > profiles, expected 


underprediction for the high lip-blowing rate case. Also shown in Fig. 36 are the root 
mean square of the optical path difference fluctuations for two additional integration 
paths. The result for the integration path which extends from r 0 = -8" to r/ = 12" 
displays an increment in < OPD' > of about (7 x 10 -7 )" from the 7" path length 
case which originates at r 0 = —3". The path initialized above the shear layer, from 
r 0 = 4" to r f = 12", shows a small < OPD' >. Finally, the time averaged optical 
path difference, OPD , can be seen to contribute curvature to the wavefronts as the 
light propagates through the shear layer. The optical clarity of the shear layer was 
determined using a (3 = 2.584 x 10 -4 , matching the value which was used to reduce 
the experimental data. 

The analytic result for the < OPD' >, which goes like x, is found from [103] 

(^^)' - 9f - ^u.,r (a* 


Derivation of the model, which utilizes time-mean quantities to determine < OPD' >, 
is given by Bogdanoff [103]. This analytic result assumes an index-matched shear layer 
with a sinusoidal n profile, n(y ) = sin ( 2L 2 J* m )- The constants, £ = 0.0091 

and L 50 „,m « 1.316* = 1.31(0.18x[m]), are empirical relations. The virtual origin of 
the shear layer is placed at x = 0" for this analysis. 
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Figure 36: 2-d treated cavity: streamwise variation of optical path difference quantities 
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Reduction of Experimental Data 


The reduction of the data obtained from experiment [104] is noted here to delineate 
the approximations used and estimate systematic error bounds in the optical path 
distortion levels. Values of p' are computed from assumptions of quasi-steady flow: 

hj = c p T ^1 -I- 


differentiation with respect to t 

Rr = PP^Pl + uu ,2^± 

P 2 7 

using {pu)' = pv! + pfu then 

S = ^ +(7 _ 1) m^-[(7-I^ + 1 ]^ 

The experimental observations against which the computed results are compared 
assume simultaneously small fluctuations in pressure and total temperature [105, 106], 
resulting in 

( 8 ) 


P V(7-1)M J pu 

Mean Mach number and density profiles are determined from isentropic relations, 
while ^=f is proportional to the voltage fluctuation, -|r, obtained from hot film probes. 
The optical path disturbance is then found from [43] 

2 rL 

( 9 ) 


{OPD'f = 2 f L <p'> 2 l r dr 

\PSTPJ Jo 


where 4- is the turbulent eddy size relative to the shear layer width, determined from 
cross correlation data to be typically about 15%. 

The few available independent measurements [106, 107] indicate that pressure 
fluctuations of about 2% of freestream static pressure occur in the shear layer spanning 
the aperture of a quieted cavity geometry. In fact, Hahn [94] reported pressure 
fluctuations of 8% from shear layer rake measurements, however these include the 
dynamic pressure component normal to the orifice as well. Pressure fluctuation levels 
can also be inferred from sound pressure levels in the cavity, observed to be at least 
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130 dB for the AOA case. Shear layer total temperature fluctuations of about 1% 
have also been reported for this Mach regime [107], The present low lip-blowing 
computation found a < p' >* 1% and a < V T >» 0.8% in the shear layer. The 
assumption of <T' r >, < p' >« 0 in a shear layer is therefore questionable, and is 
used to estimate systematic experimental error bounds. 

The determination of the error in < p' > due to background noise levels begins 
by assuming the passage of a compression wave parallel to the static pressure port in 
the wake rake. Normal reflection of the wave would impart a larger deviation from 
<_^ > as computed by Eq. 8. Utilizing Gortler’s free shear layer solution to provide 
u(r), assuming a cavity temperature recovery factor of unity, and holding mean static 
pressure constant through the layer, then p(r) is defined. The sensitivities of £ to 
| and J are ± h _^ M2+l and respectively. Using a compression or 

rarefaction wave of strength < p' > through the shear layer, then local values of p' 
due to wave passage are defined. This value of p' provides the error bound about the 
value obtained from Eq. 8, which assumes negligible < p' > and < T' t >. Taking 
shear layer pressure fluctuation levels corresponding to 135 dB and a velocity ratio 
r = 0.1, then the systematic error in the density fluctuations is 0.13% at the shear 
layer center. Figure 33 shows the resultant systematic error bars in < p' >. 

From Eq. 9 the value of < p' > is linearly proportional to < OPD' >. The error 
in < OPD' > can be found by using a conservative within-system error of 0.05% in 

pL gl eane d from Fig. 33, plus the systematic error from the above analysis. The 
resultant error bar is plotted in Fig. 36. 

5.2.5 Clean Configuration 

The SOFIA configurations were initialized using the steady solution about a clean, or 
without cavity, geometry. In order to provide a measure of validation, the geometry 
and flow conditions were chosen to replicate the wind tunnel tests: 

Moo = 0.85, Re L = 4.2 x 10 6 , L = 12.6 in. 

Poo = 0.84 kg/m 3 , p^ = 7.7 x 10 4 N/m 2 , a = 2.5° 
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However, subsequent correction of the wind tunnel data has resulted in the Ret 
4.0 x 10 6 , about 5% lower than above. The temperature difference from this correction 
causes the sound speed in these computations to be approximately 6% higher than 
experiment. Numerical results obtained for this 7% scale model are discussed below. 

The wind tunnel model without cavity was simulated in order to assess angle 
of attack errors owing to wind tunnel wall effects. Figure 37 compares the present 
pressure coefficient profiles along the crest, side, and bottom of the model with flight 
data [108] and wind tunnel results [16]. Experimental results are shown for both 
untripped and tripped cases; the latter case was used for all subsequent wind tunnel 
testing. The computations specified turbulent walls at all no-slip boundaries. Al- 
though this comparison indicates that the influence of the tunnel wall was small near 
the cavity, pressures along the bottom of the model are shifted, possibly due to the 
effect of the lower wall. A four-order drop in magnitude of 6p\ max was attained for 
this steady case in 2000 steps, using approximately four Cray Y-MP CPU hours. 


5.2.6 Configuration 25 

The geometry shown in Fig. 38 was the initial cavity configuration tested in the wind 
tunnel. This simulation was implemented in order to demonstrate the capture of 
self-excited cavity resonance in three dimensions. The flow conditions were the same 
as used above, the flowfield was initialized from the steady clean case. The stability- 
limited time step used was At = 3.53 ps. This interval size corresponds to a CFL « 1 
in the streamwise direction within the shear layer, and a C F L\ max ~ 500. 

Instantaneous Mach number contours in Fig. 39 show the flapping of the shear 
layer and interpolation treatment. Sample pressure histories on the cavity walls and 
a comparison of the PSD resulting from the wind tunnel and numerical efforts are 
shown in Fig. 40. The PSD was obtained using 2048 points, a Hanning window, and 
no zero-padding. The predicted frequencies of the dominant tones appear reasonable, 
and the computed dominant tone is within 3 dB of experiment. The magnitudes of 
the computed higher modes are much lower than observed experimentally. 

Estimation of the grid resolution required to maintain a propagating wave of a 
specific magnitude can be deduced from the rectangular two-dimensional cavity and 
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Figure 37: Clean configuration: pressure coefficient comparison 
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configuration 25 results. First, wavelength can be estimated by assuming the wave to 
be harmonic at a given frequency and travelling at the local speed of sound. Next, it 
is noted that frequencies around 2 kHz were resolved well in the two-dimensional case, 
in which the grid resolution was such that about 40 points supported the wave. From 
the configuration 25 results, it is seen that only the 700 Hz peak is well resolved, which 
again gives approximately 40 points across the wave for this coarser grid. Although 
numerical damping of the higher frequencies can be expected, most of the energy is 
contained in the lowest frequency mode, as can be seen in the sound pressure level, 



Bulkheads 

Forward Aft 



Figure 41: Comparison of sound pressure levels 
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or SPL , comparison of Fig. 41 where 

SPL MB] = 20 log,„ < l/> W™ 2 ! 

J 610 2 x 10- 5 

Experimentally observed and computed levels of < p' > in the time domain were used 
in Fig. 41. The SPL for the resonating and quieted geometries obtained numerically 
are in reasonable agreement with experiment. 

5.2.7 Configuration 100 

In 1990, an investigation of SOFIA cavity quieting treatments was performed in the 
NASA Ames 14' x 14' wind tunnel [16]. Of the many geometries tested, configuration 
100 resulted in the lowest sound production levels. This simulation was implemented 
in order to determine if the same level of quieting could be predicted numerically as 
was observed experimentally. As commented on earlier, previous investigations [41] of 
cavity noise suppression have shown aft ramp treatments to be effective by allowing a 
stable shear layer reattachment site. For the SOFIA experiment, this type of geometry 
treatment was found to be quieter than the untreated configuration 25 case by over 30 
dB. Figure 41 summarizes that the proper trends were computed. The flow conditions 
were again initialized from the clean case, and integrated using a stability-limited time 
step size of At = 7.06/xs. The frequency domain analysis was obtained using 4096 
points, a Hanning window, and no zero-padding. 

For reference purposes, Fig. 42 shows the position and orientation of the telescope 
assembly in the aircraft and the associated coarsened grids. Figure 43 shows the 
topology used in the cavity region, where the grids have again been coarsened for 
clarity. The choice of topology, driven by grid quality and turbulence modelling 
considerations, is similar to those used in the two-dimensional studies. 

Quantitative comparisons were made for this passively quieted cavity geometry in 
terms of shear layer profiles and pressure spectra in the cavity. Since errors in shear 
layer mass entrainment rate would adversely affect the cavity velocity field and hence 
the mean telescope loads, an important validation parameter is the shear layer spread 
rate. Figure 44 depicts mean experimental and computational shear layer Mach 
number profiles. The vertical scale of the profiles is twice that shown in the contour 
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Figure 42: Telescope location and grids 

plot for clarity. Note that the experimental profiles were obtained using a rake, which 
is sensitive only to u, the x-component of velocity. The discrepancy between \u\/c and 
|V|/c was found to be approximately 0.05 at the lower tail of the profile. Figure 44 
indicates reasonable agreement for growth rates, though the profile shapes become 
somewhat different as the shear layer approaches the ramp. This discrepancy may be 
in part due to probe position uncertainty and geometry modifications to allow for the 
probe mechanism. These modifications included removal of the telescope assembly 
and cutting a streamwise slot in the ramp. The difference in spread rates may also 
be due to specification of an overly-large value of o$. Time averaging of velocities 
was performed over 1000 time steps and the profiles were insensitive to the duration 
of the time-segment used. 

Some measure of qualitative agreement may be gleaned from the instantaneous 
streamlines depicted in Fig. 45, which show a strong cross flow component at the 
aft ramp for this aperture elevation angle. Although oil flow visualization was not 
performed on the configuration 100 experimental runs, a similar aft molding shown 
in Fig. 46 also displays strong cross flow behavior. 

Assessment of the oscillating telescope assembly loads requires the accurate res- 
olution of the unsteady pressure field in the cavity. Comparison of the computed 
and observed spectra at specific locations provides a measure of confidence for the 
computed telescope loads. Toward the estimation of loads, Fig. 47 shows the pres- 
sure history and resultant PSD on the cavity walls. Although the peak levels are in 
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Figure 43: Configuration 100: cavity region topology 
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Figure 44: Configuration 100: instantaneous Mach number contours and mean profiles 



Figure 45: Configuration 100: streamlines along the treated aperture 
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Figure 46: Experimentally [16] observed surface flow pattern 


agreement, the computed spectra can again be seen to drop more rapidly with fre- 
quency than the experimental results. The pressures on the primary mirror, shown 
in Fig. 48, show lower high frequency content with the magnitude of the peak at 1800 
Hz not well resolved. Figure 49 shows a low frequency component at the downstream 
secondary mirror location which was not found experimentally. The discrepancy is 
manifested as the difference between the computed and measured SPL seen in Fig. 41 
at probe 9. 

Scaling to Flight 

The above computations were performed at wind tunnel geometric scale in order to 
allow close comparison with the SOFIA experiments. An early misunderstanding 
resulted in a mismatch of the ambient conditions between the wind tunnel tests and 
the computations. Scaling of optical effects as well as the acoustic frequencies and 
magnitudes from tunnel and computation to flight are obviously of design interest, 
and hence are noted here for completeness and clarity. 









Power spectral density, PSD, dB Coefficient of Pressure, C 
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Frequency, f, Hz 


Figure 48: Configuration 100: pressure history and power spectra on primary mirror 









CHAPTER 5. RESULTS AND DISCUSSION 


101 


Briefly, the frequency is a function of the cavity length and the recovery temper- 
ature, which affects the acoustic speed in the cavity. The magnitude of the acoustic 
oscillations scale with dynamic pressure, while the optical distortion is proportional to 
density and geometric scale. Model and flight conditions are summarized in Table 2. 
Scaling from wind tunnel to flight results in a decrease in frequency by a multiplica- 



Magnitude 

Quantity 

Units 

Wind tunnel 

Computation 

Flight 


- 

0.85 

0.85 

0.85 

Rti 

- 

4 x 10 6 

4.2 x 10 6 

2.3 x 10 7 


°K 

286 

322 

217 

Poo 

±2. 

m 3 

0.77 

0.84 

0.289 

L 

m 

0.32 

0.32 

4.6 


Table 2: Wind tunnel, computed, and flight conditions 


tion factor of 16.4, a reduction of the magnitude of pressure oscillations by a factor 
^ dB, and an increase in optical distortion by 5.4 times. In contrast, scaling from 
computation to flight results in a frequency reduction of 17.4 times, a decrease in 
fluctuating pressure magnitude of 13 dB, and an increase in optical distortion of 4.9 
times. Generally, the levels of uncertainty in measured quantities is greater than the 
difference between computed and wind tunnel results. 

5.2.8 Configuration 100: Aero-Optical Effects 

The optics code was applied to the computed density field obtained for configuration 
100 from t = 0 to 7.8 ms in the manner depicted in Fig. 50. Ten rays were prop- 
agated through 110 instantaneous density fields in time intervals of At = 70. 6fis. 
Using the elapsed time for a shear layer structure to convect across the aperture as a 
characteristic time, T c = , then the optical measurement was taken for about five 

T c . The results presented here are for a computational plane at approximately the 
cross flow center of the aperture, partly shown in Fig. 50, which will provide only a 
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Station, x , inches 

Figure 50: Configuration 100: instantaneous density field and optical refraction model 

streamwise variation in optical properties. The numerical results are presented com- 
pared to previous analysis [103] and experiment [16] in which shear layer aerodynamic 
measurements were used to infer distortion. 

The levels of fluctuating density were severely underpredicted as compared to 
experiment, as can be seen in Fig. 51. Although peaks in the density fluctuations 
were computed, the highly-ordered shear layer structures similar to those found in 
the AO A study were not observed. Differences may be attributable to grid coarseness 
or within-system errors in measurements, most likely the former. 

The optical wavefront distortion through the configuration 100 aero- window is 
summarized in Fig. 52. The uppermost plot of Fig. 52 shows that the distortion model 
applied through the shear layer alone underpredicts the data determined analytically 
and experimentally. However, the computed trend is generally consistent with the 
data. At the streamwise center of the aperture, the < OPD' > at two additional 
spanwise locations are shown. These points provide an estimate of the crossflow 
variation in distortion levels. 

The center plot of Fig. 52 depicts computed < OPD' > for ray propagation 
originating below the secondary mirror, r® = —3.7", and above the shear layer, ro = 
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Figure 51: Configuration 100: density fluctuation at crossflow center of aperture 


2.3". Comparison of the computed results show an increment in < OPD' > below 
the secondary mirror. This distortion increment appears to be caused by a jet of 
re-entrant fluid originating from the shear layer impingement upon the aft ramp. 
Finally, the last plot of Fig. 52 shows that curvature is imparted to the mean optical 
field. The dip in the fluctuating and mean OPD levels at x = 42" is caused by the 
presence of the secondary mirror, in which the index of refraction, n, was fixed at 
unity. 
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Figure 52: Configuration 100: comparison of wavefront distortion 




Chapter 6 
Conclusions 


The objective of this effort was to develop and assess computational methods as ap- 
plied to unsteady perfect gas flows. The developed technology was demonstrated 
through application to two classes of problems in which unsteady effects play a dom- 
inant role. 

6.1.1 Blast- Wave Problem 

The application of two upwind schemes to unsteady, multidimensional problems 
within a structured finite-volume framework has been demonstrated on the viscous 
three-dimensional blast- wave problem. The use of time-conservative differencing and 
an approximate Riemann solver coupled with total variation diminishing methods has 
resulted in time accurate nonoscillator}' flowfield resolution. Newton subiterations are 
utilized to reduce the numerical approximations made, such as factorization error and 
the inclusion of only the first-order terms in the formation of the inviscid Jacobian. 
In addition, analysis and application of two flux evaluation methods produced only 
small differences. Finally, for the blast-wave/target interaction problem the effect of 
viscosity was increasingly significant at later times. 

Further efforts to increase the accuracy and efficiency of these methods may be 
directed along the use of nonfactored schemes or implementation on parallel ma- 
chines. Geometries of realistic complexity will require a zonal approach, necessarily 
conservative because of the strongly unsteady compressible flow regimes considered 
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here. Efficient adaptive grid techniques will reduce the memory and time expense. 
Synthesis of dynamical, structural, and fluid flow effects may provide the capability 
for an interdisciplinary simulation of the physical processes involved in this class of 
problems. 

6.1.2 Cavity Flow Problem 

The work presented here is the initial effort towards development of a cavity flow 
design and analysis tool, specifically tailored for use throughout the SOFIA project 
life. Thus far, this investigation has demonstrated that self-induced cavity resonance 
can be accurately captured for complex geometries modelled using an overset mesh 
topology. Shear layer profiles and resonant behavior are consistent with previous 
analytic and experimental work. Generally, sound pressure levels agree to within 4%. 
Topology treatment has allowed the simple specification of turbulent wall and shear 
layer regions as well as providing a means of isolating the unsteady flow region. 

Improvements in the energy distribution in frequency may be attained by use of 
higher-order spatial approximations or more simply by grid refinement. The use of 
higher order turbulence models should also be investigated. 

6.1.3 Aero-Optical Effort 

Comparison of computed and experimentally observed optical distortion levels showed 
similar trends, albeit with a discrepancy in magnitude. Large structures in the shear 
layer of the two-dimensional quieted cavity resulted in a 25% overprediction of wave- 
front distortion. The three-dimensional results underpredicted phase distortion by 
approximately 50% as compared to experiment. 

Further investigation is required to determine if improved wavefront distortion 
results can be achieved. Quantification of the effects of shear layer flow resolution on 
the optical field is certainly warranted. Improvements in the optical modelling could 
include an empirical model to account for scattering owing to subgrid turbulent scales. 
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6.1.4 Future Directions 

The development of these types of tools partially fulfills the objective of augmentation 
of experimental test programs, possibly eliminating the need to test certain specific 
configurations altogether. The use of the unsteady flowfield information appears fea- 
sible for analysis of the blast-wave/target interaction and cavity flow problem classes. 
However, current limits in computational speed prevent rapid solution throughput. 
Thus, the goal of nesting full unsteady Navier-Stokes methods within the design cycle 
awaits the advent of computers of increased power. 
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